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Abstract 

Simulation  of  reactive  energetic  materials  has  long  been  the  interest  of  researchers  due  to  its 
extensive  applications  in  military.  While  substantial  research  has  been  done  on  the  subject  at 
macro  scale,  the  research  at  micro  scale  has  initiated  recently.  Equation  of  state  (EoS)  is  the 
relation  between  physical  quantities  (pressure,  temperature,  energy  and  volume)  describing 
thermodynamic  states  of  materials  under  a  given  set  of  conditions.  It  is  the  key  connecting  the 
microscopic  and  macroscopic  behavior  when  simulating  the  macro  effects  of  an  explosion.  For 
instance,  ignition  and  growth  model  for  high  explosives  uses  two  Jones-Wilkins-Lee  (JWL) 
EoSs,  one  for  solid  explosive  and  the  other  for  gaseous  products  of  high  explosive.  The  JWL 
EoS  of  gaseous  products  is  obtained  by  conducting  explosion  experiments,  which  are  expensive, 
time-consuming  and  hazardous.  It  would  be  much  less  expensive  and  hazardous  to  obtain  the 
EoSs  of  energetic  materials  through  computational  method  instead  of  field  experiments.  In  this 
project,  we  developed  a  novel  three-stage  coupled  multi-scale  QM-MD-SPH  approach  to  model 
complex  heterogeneous  multicomponent  reactive  systems.  Our  QM-MD-SPH  approach  is 
particle-based  covering  nine  orders  of  magnitude  length  scales  (nm-m)  and  twelve  orders  of 
magnitude  timescales  (ps-s)  and  it  is  capable  of  easing  accommodation  of  multi-species  systems, 
and  each  of  the  species  as  well  as  their  interfaces  are  traceable  during  the  simulation  of  explosion 
events.  In  this  work,  the  EoSs  of  solid  /3-HMX  and  gaseous  products  of  /?-HMX  are  calculated 
using  molecular  dynamics  (MD)  simulation  with  ReaxFF-d3,  a  reactive  force  field  obtained 
from  quantum  mechanics  (QM)  simulations  based  on  the  first  principle.  The  microscopic  QM- 
MD  simulation  results  are  shown  to  compare  well  with  the  experimental  data,  other  numerical 
methods  and  some  empirical  ignition  and  growth  models.  The  calculated  QM-MD  EoSs  are  then 
incorporated  into  our  smoothed  particle  hydrodynamics  (SPH)  code  to  simulate  the  macro  effects 
of  explosions  of  high  explosive.  An  ignition  and  growth  model  is  integrated  into  the  SPH  codes 
to  simulate  the  detonation  of  high  explosives.  Further,  an  afterburning  model  is  integrated  in 
ignition  and  growth  model  and  implemented  in  our  QM-MD-SPH  package  to  simulate  explosion 
of  aluminized  explosives.  In  addition,  a  Godunov  SPH  method  is  also  incorporated  in  our  SPH 
package,  which  integrates  Riemann  solver  and  eliminates  the  need  for  artificial  viscosity  in  the 
standard  SPH.  Our  QM-MD-SPH  model  is  validated  using  benchmark  numerical  example 
problems  in  one-  and  two-dimensions  of  detonation  of  PBX  9501  explosives.  The  techniques 
proposed  in  this  report  can  also  be  used  to  simulate  the  terminal  effects  of  both  high  explosives 
(HEs)  and  casing  breakage/fragmentation. 
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Chapter  1-Introduction:  background  and  motivation 

HMX,  or  octahydro-l,3,5,7-tetranitro-l,3,5,7-tetrazocine,  is  a  powerful  high  explosive  with  extensive 
applications  in  the  military  industry.  For  example,  it  is  used  to  implode  the  fissionable  material  in  a 
nuclear  device,  as  a  component  of  plastic-bonded  explosives,  and  as  a  component  of  rocket  propellant.  It 
is  of  great  importance  to  investigate  the  properties  of  HMX  numerically  and  experimentally.  The  typical 
morphologies  of  HMX  crystals  are  shown  in  Figure  1.1  [1]. 


Figure  1.1  The  typical  morphologies  of  HMX  crystals 


The  equation  of  state  (EoS),  which  describes  the  thermodynamic  states  of  materials  under  a  given  set 
of  conditions,  plays  a  significant  role  in  determining  the  characteristics  of  energetic  materials,  such  as 
Chapman- Jouguet  (CJ)  point  and  detonation  velocity  [2].  Furthermore,  the  EoS  is  the  key  to  connect 
microscopic  and  macroscopic  phenomena  in  the  study  of  energetic  materials.  For  example,  the  EoS  is  the 
foundation  of  various  numerical  models  for  high  explosives  [3-7],  such  as  the  JWL++  model,  and  ignition 
and  growth  model.  The  JWL++  model  consists  of  two  EoSs,  one  for  the  solid  explosives  and  the  other  for 
gaseous  products  of  high  explosives.  The  two  EoSs  are  mixed  using  a  simple  additive  rule  according  to 
the  reaction  progress  of  the  explosive.  The  JWL++  model  does  not  include  any  physical  assumptions,  thus, 
it  is  only  a  primitive  numerical  model  and  the  results  produced  by  the  model  are  not  accurate.  The  other 
numerical  model,  ignition  and  growth  model  also  consists  of  two  EoSs,  similar  to  JWL++  model.  The 
difference  is  that  in  the  ignition  and  growth  model,  two  physical  assumptions  (pressure  equilibrium  and 
temperature  equilibrium)  is  introduced,  leading  to  a  more  complicated  but  also  more  accurate  mixture  rule. 
As  the  most  popular  numerical  model  for  the  simulation  of  high  explosives,  ignition  and  growth  model 
has  been  integrated  in  commercial  software  ANSYS  and  LS-DYNA. 

EoSs  of  high  explosives  are  usually  obtained  from  fitting  experimental  data.  Much  research  has  been 
done  on  the  Hugoniot  EoS  of  solid  materials.  Walsh  et  al.  initiated  the  idea  in  1955  to  use  planar  shock 
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waves  to  determine  the  EoS  of  condensed  materials  [8].  Afterwards,  shock  wave  experiment  has  been 
extensively  applied  to  obtain  the  EoS  of  a  wide  range  of  materials,  such  as  elements,  compounds,  alloys, 
porous  materials  [9-13].  Some  researchers  successfully  applied  planar  shock  waves  in  the  study  of 
Hugoniot  EoS  of  explosives  [14,  15].  Specially,  Yoo  et  al.  obtained  the  EoS  of  crystal  HMX  using 
isotropic  pressure  experiment  [16].  To  obtain  EoS  of  gaseous  products  of  high  explosives,  cylinder  test  is 
often  conducted.  Cylinder  test  involves  recording  the  detonation  velocity  and  wall  displacement  of  a 
copper  tube  filled  with  explosive  sample,  which  is  detonated  from  one  end.  The  standard  tube  is  1  in 
internal,  0.1  in  wall  and  12  in  long  [17].  The  results  of  cylinder  test  are  often  tuned  and  fitted  to  obtain 
EoS  of  gaseous  products  and  then  applied  in  numerical  models  mentioned  above. 

However,  it  is  usually  expensive  and  hazardous  to  conduct  experiments  on  high  explosives.  To  reduce 
cost  and  unnecessary  risks,  an  economic  method  is  needed  to  obtain  EoSs  of  solid  explosive  and  gaseous 
products  of  explosive.  Thus,  molecular  dynamics  (MD)  simulation  is  considered  to  fulfill  the  task. 
Molecular  dynamics  is  a  powerful  tool  to  investigate  microscopic  processes.  In  recent  years,  the 
detonation  of  high  explosives  has  been  explored  more  profoundly  with  the  introduction  of  reactive  force 
fields  (ReaxFF)  in  molecular  dynamics  [18].  Furthermore,  Liu  et  al.  introduced  an  additional  vdW-like 
interaction  to  ReaxFF,  which  enables  ReaxFF  to  obtain  the  correct  density  of  molecular  crystals  [19]. 
Strachan  et  al.  studied  thermal  decomposition  of  RDX  under  various  temperature  and  densities  [20]. 
Chenoweth  et  al.  studied  thermal  decomposition  of  a  poly  polymer  with  ReaxFF  [21].  Zhang  et  al.  studied 
thermal  decomposition  of  1,3,5-trinitrohexahydro-s-triazine  (RDX)  bonded  with  polyurethane  (Estane) 
using  molecular  dynamics  simulation  equipped  with  ReaxFF  [22].  However,  their  research  mainly 
focuses  on  the  thermal  decomposition  of  high  explosives.  In  this  work,  the  geometry  of  the  primitive  unit 
cell  of  the  P-HMX  is  obtained  from  our  previous  DFT-D2  study  [23].  The  EoSs  of  solid  explosive  and 
gaseous  products  of  p-HMX  are  calculated  using  ReaxFF-d3,  a  reactive  force  field  obtained  from 
quantum  mechanics  simulation.  The  calculated  EoSs  of  P-HMX  are  compared  with  experiment  [16]  and 
EoSs  in  ignition  and  growth  model  [24].  Then,  the  calculated  EoSs  are  incorporated  in  the  smoothed 
particle  hydrodynamics  (SPH)  code  to  simulate  the  macroscopic  explosion  of  PBX  9501,  which  contains 
95%  HMX,  2.5%  Estane  and  2.5%  BDNPA-F. 

Aluminized  explosive  is  an  enhanced  metalized  explosive  (e.g.  HMX-A1),  which  usually  contains  high 
explosive  such  as  HMX,  and  uses  aluminum  particles  (5~10  qm  in  diameter)  as  additive.  The  mass 
fraction  of  Al  powder  ranges  from  20%  to  40%  or  even  higher.  Aluminized  explosives  have  extensive 
application  in  military  industry.  Different  from  ideal  explosives  like  TNT,  HMX,  RDX,  and  TATB, 
aluminized  explosives  feature  both  fast  detonation  of  high  explosives  and  slow  combustion  of  aluminum 
particles.  The  detonation  of  high  explosives  occurs  in  microseconds,  oppositely,  the  combustion  of 
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aluminum  particles  takes  milliseconds  or  even  longer.  Due  to  the  addition  of  aluminum  powder, 
aluminized  explosives  have  relatively  low  brisance  but  high  blast  potential  [25]. 

The  detonation  of  aluminized  explosives  is  rather  complicated.  The  underlying  mechanism  of  the 
chemical  reaction  is  still  not  thoroughly  understood  till  today.  Modeling  the  detonation  and  subsequent 
combustion  requires  a  particular  knowledge  of  the  chemical  kinetics  [26]. 

To  investigate  the  detonation  of  high  explosives  and  aluminized  explosives,  researches  have 
conducted  experiments  and  performed  theoretical  analyses.  Wilkins  et  al.  obtained  the  Chapman- Jouguet 
pressure  and  equation  of  state  of  PBX  9404  through  experiments  [27].  Wackerle  et  al.  investigated  the 
pressure  history  of  PBX  9404  through  planar  shock  initiation  [28].  Cook  et  al.  measured  detonation 
pressure  and  detonation  velocity  of  TNT-A1  and  other  aluminized  explosives  [25].  Due  to  the  complexity 
of  the  problem,  it  takes  longer  time  for  researchers  to  understand  the  mechanism  theoretically.  Boiko  et  al. 
studied  the  ignition  of  aluminum  powers  in  shock  waves,  and  measured  the  relationship  between  reaction 
rate  and  temperature  [29].  Khul  et  al.  proposed  a  model  to  calculate  the  reaction  rate  of  the  combustion  of 
aluminum  under  shock  waves  [30-32].  Based  on  these  models,  Kuhl  et  al.  investigated  the  explosion  of 
aluminized  explosive  [31].  However,  the  model  proposed  by  Kuhl  et  al.  is  way  too  complicated.  It  also 
has  other  disadvantages  such  as  negative  internal  energy,  which  is  incompatible  with  mainstream  model 
of  high  explosive  such  as  ignition  and  growth  model.  Thus,  a  simplified  afterburning  model  is  employed 
in  this  work.  The  idea  is  to  use  a  reaction  rate  model  to  describe  the  combustion  of  aluminum  particles. 
Heat  will  be  gradually  released  along  with  the  chemical  reaction.  The  details  of  the  model  will  be 
introduced  in  the  following  chapters.  There  have  been  abundant  literatures  on  the  application  of 
afterburning  model  to  simulate  combustion  of  aluminum  particles.  Togashi  et  al.  investigated  the 
detonation  and  afterburning  effects  of  AFX  757  and  TNT-A1  in  confined  facility  [33,  34].  The  simulation 
agrees  well  with  measurements.  Zhou  et  al.  studied  the  detonation  of  aluminized  explosive  using  the  same 
model  and  good  agreement  with  experiment  is  observed  [35].  However,  they  just  use  simple  JWL  model 
to  simulate  the  detonation  of  high  explosives.  In  this  work,  afterburning  model  is  combined  with  ignition 
and  growth  model  so  that  the  detonation  of  high  explosive  and  combustion  of  aluminum  particles  can  be 
investigated  simultaneously. 

The  numerical  method  used  in  this  work  is  smoothed  particle  hydrodynamics  (SPH)  method.  SPH  is  a 
mesh-free  Lagrangian  method  which  was  first  developed  by  Gingold,  Monaghan  and  Lucy  et  al.  [36,  37] 
to  study  astrophysics.  Compared  with  Euler  method,  SPH  method  is  advantageous  in  tracking  free 
interfaces  and  dealing  with  large  deformation,  thus  it  has  been  extended  to  fluid  dynamics  and  solid 
mechanics  to  simulate  high-velocity  impact  and  explosion.  Liu  et  al.  [38-41]  investigated  the  feasibility  of 
using  SPH  to  simulate  explosion  of  high  explosives,  underwater  explosion,  and  shaped  charge  detonation. 
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However,  in  these  literatures,  the  numerical  model  for  explosive  is  just  simple  JWL  model,  which  is  much 
less  accurate  than  ignition  and  growth  model. 

In  this  project,  a  SPH  method  incorporated  with  ignition  and  growth  model  is  proposed  to  study  high 
explosives  and  aluminized  explosives.  The  simulation  results  are  compared  with  experiment  and  finite 
volume  method  (FVM)  solution,  and  good  agreement  is  observed,  which  indicates  that  the  proposed 
method  is  reliable  and  sufficiently  accurate.  Then,  the  EoSs  calculated  from  molecular  dynamics 
simulation  are  applied  in  SPH  codes  to  investigate  high  explosive  and  aluminized  explosive. 

A  more  advanced  SPH  method,  Godunov  method  is  described,  which  integrates  Riemann  solver  and 
eliminates  the  need  for  artificial  viscosity,  which  is  a  significant  improvement  compared  with  traditional 
SPH  with  artificial  viscosity.  A  shock  tube  problem  is  investigated  using  Godunov  SPH  to  validate  the 
accuracy  of  the  algorithm.  Then,  the  Godunov  SPH  is  applied  in  the  simulation  of  high  explosive. 

A  three-dimensional  SPH  model  with  JWL++  model  is  developed  for  non-ideal  explosive  ANFO.  The 
influence  of  smoothing  length  on  simulation  is  investigated.  A  three-dimensional  cylinder  test  model  is 
developed  to  test  the  stability  and  accuracy  of  the  proposed  SPH  method  with  JWL++  model. 
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Chapter-2  DFT-D2  study  of  HMX  using  Quantum  Mechanics  (QM)  Approach 

2.1  Accurate  modeling  of  the  density  of  HMX  using  DFT-D2  study 

Density  is  one  of  the  most  important  physical  properties  of  solid  energetic  materials  [42,  43]  because  it 
determines  the  detonation  velocity  and  detonation  pressure,  which  are  used  to  assess  the  potential 
performance  of  this  energetic  material  [44].  In  addition,  density  plays  an  important  role  in  the  equations 
of  state  describing  thermodynamic  properties,  which  can  be  used  to  describe  materials  at  continuum  level 

[45] .  Despite  its  importance,  the  accurate  prediction  of  the  densities  of  energetic  materials  is  challenging 

[46] .  Although  first-principles  calculations  based  on  density  functional  theory  (DFT)  provide  overall 
better  predictive  power  than  force  field  models,  they  fail  in  predicting  densities  of  energetic  materials 
with  standard  approximations,  partially  due  to  their  poor  descriptions  of  dispersion  forces  in  molecular 
crystals  [47-51]. 

There  are  extensive  studies  to  improve  the  modeling  of  dispersion  interactions,  or  van  der  Waals  (vdW) 
interactions.  These  studies  can  be  classified  into  two  approaches.  One  approach  is  the  modification  of  the 
exchange-correlation  functions  in  the  standard  DFT  formula  [52-54].  The  other  approach  is  to  add  an 
energy  term  for  empirical  corrections  to  the  total  energy  [55-64].  The  latter  approach  is  more  popular 
since  it  takes  much  less  computing  resources  than  the  former,  within  acceptable  accuracy.  Although  there 
are  dramatic  improvements  referring  to  standard  DFT  calculations,  previous  modellings  [49,  65]  with 
vdW  corrections  in  DFT-D1  method  [55]  still  show  considerable  errors  in  densities  compared  to 
experiments.  A  significant  improvement  in  the  prediction  of  equilibrium  densities  and  isothermal 
equations  of  state  for  hydrostatic  compressions  of  energetic  materials  at  nonzero  temperatures  was 
achieved  by  combining  van  der  Waals  corrections  with  thermal  and  zero-point  energy  corrections 
(denoted  as  DFT-D1-T  method)  [45].  However,  substantial  complexity  and  computing  demand  associated 
with  the  reported  DFT-D1-T  method  render  its  application  impractical  for  large  systems.  The  DFT-D2 
method  [57]  is  a  revised  version  of  DFT-D1  method  [55]  and  has  considerable  improvement  of  its 
predecessor  with  negligible  extra  computing  demands  comparing  to  the  standard  DFT  calculations.  Its 
application  in  studying  the  structural  properties  and  equations  of  state  of  /?- HMX  under  high  pressure 
needs  further  investigation. 

HMX,  also  known  as  cyclotetramethylene-tetranitramine,  or  l,3,5,7-Tetranitro-l,3,5,7-tetrazocane,  or 
tetrahexamine  tetranitramine,  is  an  important  secondary  explosive  that  is  most  commonly  used  in 
polymer-bonded  explosives,  and  as  a  solid  rocket  propellant.The  HMX  molecular  crystal  has  three  pure 
crystallographic  polymorphs,  as  a,  P,  and  8  phases.  The  P21/c  monoclinic  p-phase  molecular  crystal,  as 
shown  in  the  Figure  2.1(a),  is  the  thermo-dynamically  most  stable  polymorph  of  HMX  at  room 
temperature  and  has  highest  density.  A  molecule  with  a  ring-chain  structure  in  the  P-HMX  crystal  is 
shown  in  Figure  2.1(b).  Due  to  its  importance,  HMX  becomes  objective  of  extensive  studies,  both 
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experimental  and  theoretical.  It  is  generally  believed  that  the  P  phase  exists  only  below  375K  at  ambient 
pressure,  or  at  room  temperature  with  pressure  below  10  GPa. 


Figure  2.1  Geometry  of  p-HMX. 

The  unit  cell  of  |3-HMX  containing  two  HMX  molecules  with  ring-chain  structure  is  shown  in  Figure 
2.1.  The  four  carbon  atoms  (yellow)  and  four  nitrogen  atoms  form  the  ring-chain.  In  Figure  2.1(b),  a 
molecule  with  a  ring-chain  structure  is  shown.  The  molecule  is  center-symmetric. 

It  can  be  seen  from  Figure  2.1(b)  that  the  N1-N3  and  N2-N4  forms  major  and  minor  axis  of  the  ring- 
chain  structure,  respectively.  The  four  carbon  atoms  are  coplanar  on  theC4-plane.  The  four  nitrogen 
atoms  on  the  ring-chain  are  also  coplanar  on  the  N4-plane.The  angle  between  C4-plane  and  N4-plane  is 
30.3°. 

Here  we  model  P-HMX  using  DFT-D2  method  [57],  which  describes  the  van  der  Waals  interactions  as 
a  simple  pair-wise  force  field.  This  method  is  chosen  as  a  result  of  compromise  between  two  opposite 
considerations:  accuracy  and  feasibility.  We  investigate  the  molecular  structure,  mechanical  properties, 
and  equations  of  state  of  the  P-HMX.  To  the  best  of  our  knowledge,  for  the  first  time,  we  predict  the 
elastic  constants  of  P-HMX  using  DFT-D2  calculations. 

We  consider  a  conventional  unit  cell  containing  two  HMX  molecules  (56  atoms  in  total)  with  periodic 
boundary  conditions,  as  depicted  in  Figure  2.1.  The  total  energies  of  the  system,  forces  on  each  atom, 
stresses,  and  stress-strain  relationships  of  p-HMX  under  the  desired  deformation  configurations  are 
characterized  via  first-principles  calculations  based  on  density-functional  theory  (DFT).  DFT  calculations 
were  carried  out  with  the  Vienna  Ab-initio  Simulation  Package  (VASP)  [66,  67]  which  is  based  on  the 
Kohn-Sham  Density  Functional  Theory  (KS-DFT)  [68,  69]  with  the  generalized  gradient  approximations 
as  parameterized  by  Perdew,  Burke,  and  Ernzerhof  (PBE)  for  exchange-correlation  functions  [70].  The 
electrons  explicitly  included  in  the  calculations  are  the  Is1  for  hydrogen  atoms,  2s22p4  electrons  for 
carbon  atoms,  2s22p5  for  nitrogen  atoms,  and  2s22p6  for  oxygen  atoms.  The  core  electrons  are  replaced  by 
the  projector  augmented  wave  (PAW)  and  pseudo-potential  approach  [71,  72].  The  kinetic-energy  cutoff 
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for  the  plane-wave  basis  was  selected  to  be  800  eV  in  this  study.  The  calculations  are  performed  at  zero 
temperature. 

The  criterion  to  stop  the  relaxation  of  the  electronic  degrees  of  freedom  is  set  by  total  energy  change  to 
be  smaller  than  0.000001  eV.  The  optimized  atomic  geometry  was  achieved  through  minimizing 
Hellmann-Feynman  forces  acting  on  each  atom  until  the  maximum  forces  on  the  ions  were  smaller  than 
0.001  eV/A.  The  atomic  structures  of  all  the  deformed  and  undeformed  configurations  were  obtained  by 
fully  relaxing  a  168-atom-unit  cell.  The  simulation  invokes  periodic  boundary  conditions  for  the  three 
directions.  The  irreducible  Brillouin  Zone  was  sampled  with  a  5  X  3  X  4  Gamma-centered  k- mesh.  The 
initial  charge  densities  were  taken  as  a  superposition  of  atomic  charge  densities. 

In  the  DFT-D2  method  [57],  the  van  der  Waals  interactions  are  described  using  a  pair-wise  force  field. 
Such  a  semi-empirical  dispersion  potential  is  then  added  to  the  conventional  Kohn-Sham  DFT  energy  as 
Edft-D2  =  Edft  +  Edisp  >  and 

Edisp  =  Eqn2.1 

where  N  is  the  number  of  atoms.  The  summations  go  over  all  atoms  and  all  translations  of  the  unit  cell 
—  h>  h)  •  The  prime  indicates  that  for  =  0  ,  i  ^  j  .  C6ij-  stands  for  the  dispersion  coefficient  for  the 
atom  pair  j  .  r^L  is  the  distance  between  atom  i  in  the  reference  cell  L  —  0  and  atom  j  in  the  cell  L. 
f(jij )  is  a  damping  function  whose  role  is  to  scale  the  force  field  such  as  to  minimize  contributions  from 
interactions  within  typical  bonding  distances  r.  Since  the  van  der  Waals  interactions  decay  quickly  in  the 
power  of  -6,  the  contributions  outside  a  certain  suitably  chosen  cutoff  radius  are  negligible.  The  cutoff 
radius  for  pair  interaction  in  this  study  is  set  to  30.0  A.  Here  Fermi-type  damping  function  is  used  as 

fd,6  {jij )  —  /  rtj  T  Eqn  2.2 

where  S6  is  the  global  scaling  parameter.  The  global  scaling  factor  S6  =  0.75  is  used  for  PBE  exchange- 
correlation  functions.  SR  is  fixed  at  1.00.  The  damping  parameter  d  =  20.0  is  used. 

Compared  to  DFT-D2  method,  the  DFT-D1  method  [55]  has  a  few  shortcomings:  (1)  It  has  parameters 
available  only  for  elements  H,  C-Ne.  (2)  There  are  systematic  errors  for  calculations  of  molecules 
containing  third-row  elements.  (3)  It  leads  to  inconsistencies  for  thermochemistry  due  to  double-counting 
problem  [57].  In  addition  to  overcome  these  shortcomings,  DFT-D2  generalized  the  dispersion  correction 
and  improved  the  accuracy. 

We  first  optimize  the  geometry  of  the  monoclinic  crystal  (also  shown  in  Figure  2.1)  by  full  relaxation 
of  all  the  atoms  and  lattice  constants.  The  optimized  lattice  constants  are:  a-6.542  A,  b- 10.842  A,  c- 
8.745  A,  a=y= 90°,  and  /?=124.413°.  For  the  comparison,  we  also  optimize  the  structure  without  the  van 
der  Waals  corrections.  Our  results  of  the  lattice  constants,  the  volume  of  the  unit  cell,  and  the  densities  are 
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summarized  in  Table  2.1  and  compared  to  the  experiment  and  previous  theoretical  results.  It  shows  that 
the  standard  DFT  calculations  give  poor  predictions.  For  example,  the  volume  of  the  unit  cell  is  6.7% 
higher  than  the  experimental  value.  The  prediction  of  the  density  of  /?  -HMX  from  our  DFT-D2 
calculations  is  more  accurate,  with  a  difference  of  -1.47%  compared  to  the  experimental  value.  This  is  a 
significant  improvement  over  standard  DFT  calculations  without  van  der  Waals  corrections. 


Table  2.1  Lattice  constants  a,  b,  c,  lattice  angle  /?  ,  volume  of  the  unit  cell  V  ,  and  density  p  predicted 
from  DFT  and  DFT-D2  calculations,  compared  with  experiments  and  previous  calculations.  The  numbers 
in  parentheses  are  differences  in  percentage  referring  to  the  experiment. 


a  (A) 

b  (A) 

c  (A) 

A 

V(A3 ) 

p  (103  Kg /m3) 

Expt.a 

6.54 

11.05 

8.70 

124.30 

519.39 

1.894 

Expt.b 

6.537 

11.054 

8.7018 

124.443 

518.558 

1.897 

Expt.c 

6.5255 

11.0369 

7.3640 

102.670 

517.45 

1.901 

Expt.d 

6.54 

11.05 

7.37 

102.8 

519.37 

1.894 

DFT 

6.673(+2.03%) 

11.312(+2.37%) 

8.894(+2.23%) 

124.395(+0.08%) 

553.99(+6.66% 

1.775(-6.23%) 

DFT-D2 

6.542(+0.03%) 

10.842(-1.88%) 

8.745(+0.52%) 

124.41(+0.09%) 

511.73(-1.47%) 

1.923(+1.53%) 

Theory  e 

6.70(+2.45%) 

11.35(+2.71%) 

8.91  (+2.41  %) 

124.13(-0.14%) 

560.86(+7.98%) 

1.754(-7.39%) 

Theory  fl 

6.38(-2.45%) 

10.41(-5.79%) 

8.43(-3.1%) 

123.0(-1.05%) 

463.1(-10.74%) 

2.122(12.03%) 

Theory  f2 

6.90(+5.50%) 

11.65(+5.43%) 

9.15(+5.17%) 

124.5(+0.16%) 

608.1(+17.08%) 

1.617(-14.59%) 

Theory  f3 

6.56(+0.31%) 

10.97(-0.72%) 

8.70(0.0%) 

124.4(+0.08%) 

517.4(-0.38%) 

1.901(+0.38%) 

Theory  g2 

6.78(+3.67%) 

11.48(+3.89%) 

9.19(+5.63%) 

125.02(+0.58%) 

585.57(+12.74%) 

1 .680(-l  1 .30%) 

Theory  gl 

6.43(-1.68%) 

10.34(+6.43%) 

8.61(-1.03%) 

124.23(-0.06%) 

473.81(-8.78%) 

2.076(+9.62%) 

Theory  h 

6.539(-0.02%) 

11.03(-0.18%) 

8.689(-0.13%) 

123.9(-0.32%) 

520.17(+0.15%) 

1.891(-0. 15%) 

Theory  i 

6.762(+3.39%) 

11.461(+3.72%) 

8.865(+1.90%) 

123.8(-0.40%) 

570.599(+9.89%) 

1.724(-8.97%) 

Theory  j 

6.67(+1.99%) 

11.17(+1.09%) 

8.95(+2.87%) 

124.5(+0.16%) 

549.30(+5.76%) 

1.791(-5.45%) 

Theory  k 

6.57(+0.46%) 

11.02(-0.27%) 

9.04(+3.91%) 

124.9(+0.48%) 

531.12(+2.26%) 

1.852(-2.21%) 

Theory  1 

6.58(+0.61%) 

11.12(+0.63%) 

8.76(+0.69%) 

124.3(+0.0%) 

529.8(+2.0%) 

1.856(-1.96%) 

Theory  m 

6.57(+0.46%) 

10.63(-3.80%) 

9.13(+4.94%) 

123.67(-0.51%) 

530.6(+2.16%) 

1.854(-2.11%) 

Theory  nl 

- 

- 

- 

- 

556.07(+7.06%) 

1.769(-6.60%) 

Theory  n2 

- 

- 

- 

- 

500.77(-3.58%) 

1.964(+3.72%) 

Theory  n3 

- 

- 

- 

- 

519.41(+0.0%) 

1.894(0.00%) 

a:  Ref.73,16;  b:  Ref.74;  c:  303  K  space  group  P  21/n  in  Ref.50;  d  :  space  group  P21/n  in  Ref.75;  e:  DFT  study  using  PAW-PBE 
(GGA)  in  Ref. 49;  ft:  DFT  study  using  USPP-LDA  in  Ref.76;  f2:  DFT  study  using  USPP-PBE(GGA)  in  Ref.76;  £3:  DFT-D2 


study  using  USPP-PBE(GGA)  in  Ref.76;  gl:  DFT  study  using  USPP-LDA  in  Ref.77;  g2:  DFT  study  using  USPP-PBE(GGA)  in 
Ref.77;  h:  DFT  study  using  USPP-LDA  in  Ref.78;  i:  DFT  study  using  USPP-PBE(GGA)  in  Ref.79;  j:  Monte  Carlo  calculations 
in  Ref.80;  k:  Molecular  Dynamics  study  in  Ref.81;  1:  Molecular  Dynamics  study  in  Ref.82. 
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The  one  non-standard  cell  commonly  used  is  P2i/n  which  is  derived  from  P2i/C.  The  lattice  vectors  a0  , 
b0  ,  and  c0  in  space  group  P2i/n  can  be  transformed  to  space  group  P2i/c  as  a  =  —a0  ,  b  =  —b0  ,  and 
c  =  a0  +  b0  . 

In  general,  the  molecular  dynamics  (MD)  simulations  predict  better  lattice  vectors  than  the  standard 
DFT  calculations.  It  is  partially  because  the  MD  calculations  include  the  van  der  Waals  interactions.  Once 
the  van  der  Waals  corrections  are  introduced,  the  first-principles  calculations  [45]  show  good  predictions 
with  accuracy  within  3%  of  the  experimental  values. 


2.  2  Equation  of  State  of  HMX  up  to  100  GPa 


(a)  The  EoS  of  the  solid  P-HMX  at  zero  temperature. 


(b)  The  detailed  EoS  of  the  solid  p-HMX  at  zero  temperature.  Volume  ranges  from  400  A3  to  560  A3 


Figure  2.2  Equation  of  state.  The  pres  sure- volume  relationship  of  the  solid  P-HMX  at  zero 

temperature. 


18 


We  study  the  isothermal  equations  of  states  of  /?-HMX  at  zero  temperature  under  hydrostatic  pressures. 
The  corresponding  volume  is  obtained  after  the  system  is  fully  relaxed.  The  pressure-volume  curve  of 
unreacted  /?- HMX  at  the  temperature  of  OK  was  illustrated  in  Figure  2.2.  The  volume  corresponds  to  the 
56-atom-unitcell  volume.  The  isothermal  hydrostatic  equation  of  state  of  (3  -HMX  is  compared  with 
experiments  (Gump  [83]  and  Yoo  [16])  and  previous  calculations  (Conroy  [49]  and  Landerville  [45])  . 
The  upper  panel  shows  pressure  from  0-100  GPa.  The  corresponding  volume  in  the  present  DFT-D2 
calculations  varies  from  512.64  A3  (100%)  to  265.95  A3  (51.9%).  The  lower  panel  shows  pressure  from 
0-10  GPa  for  better  comparison.  Unlike  standard  DFT  calculation  (blue-circle  line),  our  DFT-D2  study 
(red-square  line)  shows  reasonably  good  agreement  with  the  hydrostatic-compression  experiments  [83, 
16],  suggesting  that  the  Van  der  Waals  interaction  is  critically  important  in  modeling  the  mechanical 
properties  of  this  molecular  crystal. 

In  Figure  2.2,  the  volume  is  the  56-atom-unit-cell  volume.  The  upper  panel  shows  pressure  from  0-100 
GPa.  The  corresponding  volume  in  the  present  DFT-D2  calculations  varies  from  512.64  A3  (100%)  to 
265.95  A3  (51.9%).  The  lower  panel  shows  pressure  from  0-10  GPa  for  better  comparison.  The 
experiments  (Gump  [84]  and  Yoo  [16])  and  previous  calculations  (Conroy  [49]  and  Landerville  [45]). 

It  is  worthy  pointing  out  that  the  experimental  data  were  collected  at  room  temperature.  Whereas,  our 
results  are  for  zero  temperature  and  we  have  not  corrected  the  results  to  account  for  finite  temperature. 
Moreover,  our  calculations  are  performed  using  the  (3  polymorph  of  HMX,  which  is  consistent  with  the 
experimental  data  of  Gump  and  Peiris  [83].  As  observed  in  previous  studies  [49,  84]  the  calculated 
isotherm  appears  to  approach  experimental  curve  with  increasing  pressure.  However,  as  discussed  below, 
the  sample  from  the  experiment  of  Yoo  and  Cynn  [16]  is  no  longer  in  the  (3  phase  for  pressures  beyond  12 
GPa. 

In  order  to  determine  the  bulk  modulus  and  its  derivative  with  respect  to  pressure,  the  calculated 
hydrostatic-compression  data,  up  to  about  12  GPa,  were  fit  to  the  third-order  Birch-Murnaphan  (BM) 
equations  of  states,  [85,  49]  as 

P  =  3-B0(t]-7/3  -  7 r5/3)  [l  +  l  (B'o  -  4)0T2/3  -  1)]  Eqn.  2.3 

where  77  =  V /V0  is  the  ratio  of  volume  V  in  response  to  P  to  the  equilibrium  volume  V0  at  zero  pressure. 

dB 

B0  is  the  bulk  modulus  and  B'0  =  is  the  pressure  derivative  of  bulk  modulus. 

The  values  of  B0  and  50  from  both  standard  DFT  and  DFT-D2  calculations  are  listed  in  Table  2.2, 

compared  to  experimental  values  [16,  83]  and  theoretical  predictions  [49,  45]  that  are  also  obtained  from 
fitting  to  the  BM  equations  of  states.  The  errors  from  the  fitting  process  of  this  study  are  also  listed  in  the 
Table  2.2.  Our  results  agree  with  experiments  and  the  previous  predictions.  Moreover,  our  bulk  modulus 
from  this  BM  EOS  fitting  also  agrees  with  elastic  constants  calculations  in  previous  subsection. 
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Table  2.2  The  Bulk  modulus  B  and  its  derivative  with  respect  to  pressure  B'  fitted  from  the  isothermal 
equations  of  states  from  DFT  and  DFT-D2  studies,  compared  to  experiments  and  theoretical  predictions. 

The  error-bar  rises  in  the  least-square  fitting  process. 


Source 

B0  (GPa) 

error 

Bo 

error 

Ref.  [16] 

Expt 

16.7 

6.8 

Ref.  [83] 

Expt 

21.0 

7.45 

DFT-D2 

Theory 

14.46 

±0.98 

9.96 

±0.98 

DFT 

Theory 

13.62 

±2.8 

6.5 

±0.8 

Ref.  [49] 

Theory 

9.8 

9.1 

Ref.  [45] 

Theory 

16.74 

7.78 

By  comparing  the  predictions  from  DFT-D2  and  standard  DFT,  one  can  conclude  that  the  van  der 
Waals  correction  is  necessary  for  accurate  modeling  the  elastic  properties  of  P-HMX. 

2.3  Structural  response  of  HMX  to  hydrostatic  compression:  molecular  “hot  spots” 

We  studied  the  evolution  of  the  structures  of  P-HMX  under  high  pressures.  Firstly  we  studied  the 
lattice  structures  including  lattice  constants  and  lattice  angles.  We  observed  that  the  lattice  constant  b  and 
lattice  angle  P  are  more  sensitive  to  the  applied  pressures.  We  then  studied  the  bond  lengths  of  C-H,  C-N, 
N-O,  and  N-N  bonds.  We  noticed  that  the  N-N  bonds  are  vulnerable  to  pressure.  We  further  studied  the 
bending  angle  between  N-N  bonds  with  the  plane  formed  by  carbon  atoms. We  found  the  N-N  bonds 
along  the  minor  axis  of  the  ring  are  more  susceptible  to  pressure. 

The  lattice  constants  a ,  b ,  and  c  as  a  function  of  hydrostatic  pressure  p  are  plotted  in  Figure  2.3, 
compared  with  experimental  results  from  Gump  [83]  and  Yoo  [16].  For  a  good  comparison,  three  panels 
for  the  range  from  (a)  0-10  GPa,  (b)  0-40  GPa,  and  (c)  0-100  GPa  are  presented,  while  keeping  the  y  in 
the  same  scale. There  is  a  good  agreement  of  the  lattice  constants  between  our  DFT-D2  predictions  with 
the  experiments  under  the  pressure  of  10  GPa  where  there  is  no  phase  transition.At  the  pressure  of  12 
GPa,  there  is  a  phase  transition  in  experiment. 

The  lattice  constants  a,  b  and  c  as  a  function  of  hydrostatic  pressure  p  ranging  from  0-10  GPa,  0-40 
GPa,  and  0-100  GPa  is  shown  in  Figure  2.3(a),  (b),  (c)  respectively.  Simulation  results  are  compared 
with  experiments  from  Gump  [83]  and  Yoo  [16]. 

At  the  pressure  of  36  GPa,  there  is  another  phase  transition  observed  in  the  experiment.  However, 
those  phase  transitions  were  not  observed  in  our  computations,  partly  due  to  the  small  cells  which  exclude 
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the  long  wavelength  phonon  modes  in  addition  to  the  zero  temperature  that  preclude  the  dynamic 
instabilities. 

We  noticed  that  the  lattice  constant  b  is  more  sensitive  to  the  pressure  than  lattice  constants  a  and  c. 


(a)  Lattice  constants  a,  b  and  c  as  a  function  of  hydrostatic  pressure  p  ranging  from  0-10  GPa 


Pressure  (GPa) 


(b)  Lattice  constants  a,  b  and  c  as  a  function  of  hydrostatic  pressure  p  ranging  from  0-40  GPa 
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(c)  Lattice  constants  a,  b  and  c  as  a  function  of  hydrostatic  pressure  p  ranging  from  0-100  GPa 
Figure  2.3  Lattice  constants  a,  b ,  and  c  under  high  pressure. 

Similar  to  the  lattice  constants,  we  studied  the  lattice  angles  as  a  function  of  hydrostatic  pressure.  We 
found  that  the  angle  a  and  y  varies  little  under  pressures.  On  the  contrary,  the  angle  (3  is  more  sensitive  to 
pressure.  The  (3  angle  as  a  function  of  hydrostatic  pressure  p  are  plotted  in  Figure  2.4,  and  compared  with 
experimental  results  from  Gump  [83]  and  Yoo  [16]. 

Notice  that  the  phase  transition  is  clearly  reflected  in  the  experimental  plots,  whereas,  it  is  absent  in 
the  simulation  plots  due  to  the  small  cell  and  zero  temperature  constrain  in  the  DFT-D2  modeling. In  spite 
of  these  differences,  there  is  a  reasonably  good  agreement  between  the  lattice  constants  obtained  in  our 
DFT-D2  calculations  and  experiments.  Furthermore,  we  noticed  that  the  lattice  angle  (3  is  the  most 
sensitive  to  pressure  among  all  three  lattice  angles. 

The  lattice  angle  [3  as  a  function  of  hydrostatic  pressure  p  ranging  from  0-10  GPa,  0-40  GPa,  0-100 
GPa  is  shown  in  Figure  2.4(a),  (b),  (c)  respectively.  Simulation  results  are  compared  with  experiments 
from  Gump  [83]  and  Yoo  [16]. 
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(a)  Lattice  angle  (3  as  a  function  of  hydrostatic  pressure  p  ranging  from  0-10  GPa 
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(b)  lattice  angle  (3  as  a  function  of  hydrostatic  pressure  p  ranging  from  0-40  GPa 
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(c)  lattice  angle  (3  as  a  function  of  hydrostatic  pressure  p  ranging  from  0-100  GPa 
Figure  2.4  Lattice  angle  /?  under  high  pressure. 


Next,  we  examine  the  bond  lengths  under  various  pressures  in  order  to  find  the  atomistic  mechanism 
corresponding  to  the  variations  of  lattice  structures. The  bond  lengths  of  the  bonds  C-H,  C-N,  N-O,  and  N- 
N  as  a  function  of  hydrostatic  pressure  p  ranged  from  0-100  GPa  are  illustrated  in  Figure  2.5.  In  our  unit 
cell,  the  number  of  bonds  is  16,  16,  16,  and  8  for  C-H,  C-N,  N-O,  and  N-N,  respectively.  The  bond 
lengths  are  averaged  over  the  unit  cell.  Our  results  show  that  the  N-N  bonds  are  the  most  sensitive  to  the 
applied  hydrostatic  pressure,  which  indicates  that  it  is  vulnerable  to  compression. 

The  bond  lengths  of  C-H,  C-N,  N-O,  and  N-N  as  a  function  of  hydrostatic  pressure  p  ranging  from  0- 
100  GPa  are  shown  in  Figure  2.5.  The  bond  lengths  are  averaged  over  the  unit  cell. 

In  order  to  find  the  atomic  mechanism  corresponding  to  the  variations  of  lattice  structures,  we  further 
study  examine  the  bending  angles  under  various  pressures.  There  are  two  kinds  of  N-N  bonds  in  a  HMX 
molecule:  one  along  the  major  axis  along  N1-N3  of  the  ring-chain  (Figure  2.1b)  and  the  other  along  the 
minor  axis  along  N2-N4  of  the  ring-chain  (Figure  2.1b).  The  four  carbon  atoms  are  co-planar,  marked  as 
C4-plane.  Due  to  the  symmetry,  there  are  two  angles  between  the  N-N  bonds  and  the  C4-plane.The  angle 
along  the  major  axis  is  denoted  as  pi  and  the  angle  along  the  minor  axis  is  denoted  as  P2.  During  the 
compression,  the  two  angles  piand  P2might  respond  differently,  causing  the  anisotropy  of  the  crystal. 
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Figure  2.6  Bond  angles  under  high  pressure. 

The  angles  (pi  and  P2)  between  the  N-N  bond  and  the  plane  formed  by  the  four  carbon  atoms  as  a 
function  of  hydrostatic  pressure  p  ranging  from  0-100  GPa  is  shown  in  Figure  2.6.  We  found  that  pi 
increases  with  respect  to  an  increasing  pressure,  opposite  to  the  decrease  of  [32.  Furthermore,  the  angle  P2 
is  more  sensitive  to  the  applied  pressure  than  the  angle  pi,  indicating  that  the  N-N  bonds  along  the  minor 
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axis  are  more  vulnerable  to  the  compression.Therefore,  we  may  conclude  that  the  N-N  bonds  along  the 
minor  axis  are  responsible  for  the  sensitivity  of  P-HMX. 

2.4  Elastic  properties  of  HMX 

We  next  calculate  the  elastic  properties  of  p-HMX.  The  elastic  constants  are  the  second  derivatives  of 
the  total  energy  with  regard  to  strain.  For  the  first  time,  we  predict  the  elastic  constants  of  P-HMX  using 
DFT-D2  method.  Our  results  of  all  13  elastic  constants,  in  standard  Voigt  notation  45,  are  listed  in  Table 
2.3,  and  compared  with  experiments  and  theoretical  pre-  dictions.  Our  results  show  that  the  prediction 
from  standard  DFT  studies  is  far  off  the  experiments,  but  the  prediction  from  DFT-D2  agrees  reasonably 
well  with  the  experiments. 


Table  2.3  13  Non-zero  elastic  constants  Qy,  bulk  modulus  B  ,  shear  modulus  G  ,  Poisson’s  ratio  v  , 


and  Cauchy  pressure  CP(CP  =  C12  —  C44)  of  /?- HMX  predicted  from  DFT  and  DFT-D2  calculations, 


compared  with  experiments  and  previous  calculations.  The  units  are  GPa  except  Poisson  ratio 


Ci i 

^12 

^13 

^15 

C22 

C23 

^25 

C33 

C35 

C44 

C\6 

C55 

Ce  6 

B 

G 

V 

Expt  a. 

20.58 

9.65 

9.75 

-0.61 

19.69 

12.93 

4.89 

18.24 

1.57 

9.92 

4.42 

7.69 

10.67 

13.69 

7.40 

0.271 

-0.27 

Expt  b. 

18.41 

6.37 

10.50 

-1.10 

14.41 

6.42 

0.83 

12.44 

1.08 

4.77 

2.75 

4.77 

4.46 

10.20 

4.26 

0.317 

1.6 

Expt  c. 

19.4 

5.9 

8.4 

-1.1 

17.5 

8.2 

3.2 

17.8 

0.2 

9.1 

2.4 

9.2 

9.8 

11.08 

7.77 

0.217 

-3.2 

Expt  d. 

20.8 

4.8 

12.5 

-0.5 

26.9 

5.8 

-1.9 

18.5 

1.9 

4.2 

2.9 

6.1 

2.5 

12.49 

5.43 

0.310 

0.6 

DFT 

75.1 

18.6 

1.4 

-23.4 

28.9 

24.6 

1.0 

44.5 

-14.2 

23.4 

-8.0 

18.7 

28.7 

26.41 

21.09 

0.185 

-4.8 

DFT-D2 

29.3 

10.6 

13.8 

-2.1 

25.0 

16.6 

6.2 

27.5 

1.1 

13.6 

6.8 

12.8 

13.9 

18.20 

10.78 

0.253 

-3.0 

MDf 

22.2 

9.6 

13.2 

-0.1 

23.9 

13.0 

4.7 

23.4 

1.6 

9.2 

2.5 

11.1 

10.1 

15.68 

8.33 

0.274 

0.4 

MD  g 

13.6 

3.75 

4.66 

-0.15 

9.50 

5.07 

-2.71 

13.2 

-0.96 

6.41 

-2.1 

4.04 

4.68 

7.03 

4.55 

0.234 

-2.7 

DFT  h 

22.2 

5.0 

11.5 

1.3 

23.1 

8.2 

2.2 

15.3 

3.5 

4.1 

4.7 

5.5 

1.6 

12.22 

4.63 

0.332 

0.9 

a:  Ref[75];  b:  Ref[86];  c:  Ref[87];  d:  Ref[88];e:  Ref[89];f:  Ref[82],g:  Ref[90];h:  Ref[91] 


For  a  valid  comparison,  the  Cartesian  system  to  which  the  elastic  constants  referenced  must  be  the 
same  for  all  these  data  sets.  We  used  the  reference  system  of  the  x  and  y  axes  parallel  to  the  a  -  and  b  - 
crystallographic  axes,  respectively,  which  is  the  same  as  Sewell  Bibliography  entry  and  Stevens 
Bibliography  entry  .  Zaug’s  results  ( d  in  Table  2.3)  were  transformed  to  this  Cartesian  system  by  Sewell. 
Note  that  the  experimental  studies  prefer  using  the  P21  /  n  space  group  over  standard  P21  /  c  space 
group.  The  elastic  constants  in  these  two  space  groups  are  identical  except  that  C45  =  —C4S  ,  C25  =  ~C25  , 
C25  =  —C25  ,  and  C46  =  — C46  ,  where  supper  index  c  refers  to  space  group  P21  /  c  and  super  index  n 
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refers  to  space  group  P21  /  n  .  For  a  valid  comparison,  we  transformed  our  results  of  C-  into  the  form  of 
C/}  (Table  2.3). 

The  bulk  modulus  B  measures  the  resistance  to  hydrostatic  compression  and  it  can  be  calculated 
using  the  Voigt  notation  from  the  elastic  constants  as: 

B  =  -9Z.ljCij  Eqn  2.4 

Our  calculated  value  of  the  bulk  modulus  is  18.2  GPa,  which  is  larger  than  the  average  value  11.9  GPa 
of  the  four  sets  of  experimental  values.  It  is  well  known  that  the  theoretical  elastic  constants  are  often 
stiffer  compared  to  those  of  experiments.  This  overestimation  is  attributed  to  two  factors:  the  use  of  the 
rigid-body  approximation  for  flexible  molecules,  and  the  neglect  of  the  anharmonic  softening  of  the 
lattice  at  finite  temperatures.  Our  DFT  studies  are  performed  at  zero  temperature  where  the 
thermodynamics  as  well  as  the  anharmonic  softening  are  excluded.  This  is  reflected  in  the  larger  value  of 
our  predicted  bulk  modulus. 

The  shear  modulus  measures  the  resistance  in  shearing.  The  shear  deformations  have  been  proposed  to 
be  a  critical  component  in  the  initiation  of  detonation.  Therefore,  a  precise  determination  of  the  shear 

modulus  is  necessary  for  substantiating  possible  detonation  mechanisms.  The  shear  modulus  G  can  be 
calculated  using  the  Voigt  notation  from  elastic  constants  as: 

G  =  ~  (Cn  +  C22  +  C33  -  C12  ~  C23  -  C13)  +  -  (C44  +  C55  +  C66)  Eqn  2.5 

Our  result  of  the  shear  modulus  is  10.8  GPa,  which  is  larger  than  the  average  value,  6.2  GPa,  of  the 
four  sets  of  experimental  values.  This  is  partially  attributed  to  the  preclusion  of  the  anharmonic  softening 
of  the  lattice  with  temperature. 

The  Poisson  ratio  measures  the  shape  change  under  mechanical  loading.  For  an  isotropic  material,  the 


Poisson  ratio  V  is  related  with  bulk  modules  and  shear  modulus  as: 

3B-2G  _ 

v  =  — - -  Eqn  2.6 

2(3  B  +  G)  H 

We  predicted  that  the  Poisson  ratio  is  0.253,  which  is  slightly  smaller  than  the  average  (0.278)  of  the 
four  sets  of  experimental  values.  Interestingly,  opposite  to  the  bulk  modulus  and  shear  modulus,  the 
Poisson  ratio  agrees  well  with  the  experiments. 

The  Cauchy  pressure  CP  =  C12  —  C4  4  can  be  used  to  evaluate  the  deductibility  of  a  material  [92].  The 
positive  value  of  CP  indicates  that  the  material  is  ductile  whereas  the  negative  value  implies  a  brittle 
material.  The  larger  positive  value  means  that  the  material  is  more  ductile.  The  results  of  ductibility  from 
both  experiments  and  molecular  dynamics  studies  are  controversial.  For  example,  two  experiments  [75, 
87]  indicate  /?-HMX  is  brittle,  but  other  two  [86,  88]  imply  that  /?- HMX  is  ductile.  The  average  Cp  of 
the  four  experiments  is  -0.32.  Our  DFT  and  DFT-D2  studies  indicate  that  the  /?- HMX  is  brittle,  but  with 
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much  larger  magnitudes.  The  possible  reasons  of  that  our  model  predicts  much  lower  ductibility  than 
experiments  are  twofold.  First,  the  simulation  boxes  in  the  DFT  and  DFT-D2  modeling  are  small,  which 
exclude  the  long  wavelength  phonon  modes  that  contribute  to  the  ductility.  Second,  the  zero  temperature 
applied  in  our  models  precludes  the  thermodynamics  effects,  which  is  essential  for  the  ductility.  As  a 
result,  the  ductibility  of  our  models  appears  to  be  much  smaller. 

The  diagonal  elements  Cit  describe  the  crystal  stiffness  under  uniaxial  compression  (i  =  1,2,3)  and 
shear  ( i  =  4,5,6),  whereas  the  off-diagonal  elements  C^j  presents  the  crystal  stiffness  under  biaxial 
compression  (i  ^  j  =  1,2,3)  and  distortion  (i  ^  j  =  4,5,6).  Our  results  show  that  there  is  a  remarkable 
anisotropy  in  the  diagonal  elements  of  the  elastic  constant  tensor.  Furthermore,  it  indicates  that  the/?- 
HMX  is  anisotropic  upon  compression  and  has  a  good  stability  to  shear  deformation  perpendicular  to  the 
c-axis. 

Here  C22  is  relevant  to  the  “doorway  mode”  anharmonic  coupling  model,  which  is  an  important 
feature  of  detonation  of  secondary  explosives.  [93]  Experimental  and  theoretical  vibrational  analysis 
illustrates  the  presence  of  low-energy  modes  that  have  the  necessary  symmetry  for  such  anharmonic 
coupling.  The  experimental  value  of  C22  varies  from  14.41  GPa  to  26.9  GPa,  with  an  average  of  19.63 
GPa.  Our  predicted  value  (25  GPa)  of  C22  agrees  with  experiments.  Especially,  it  is  very  close  to  the 
experimental  findings  (26.9  GPa)  [88]. 

2.5  Bandgaps  of  HMX  under  hydrostatic  compression 

The  impact  sensitivity  is  a  very  important  characteristic  of  an  energetic  material,  and  it  is  closely 
related  to  its  production,  storage,  transportation,  and  detonation.  It  is  believed  that  the  electronic  bandgap 
is  correlated  with  the  impact  sensitivity  of  an  energetic  material  [94,  78].  The  bandgap  is  the  electronic 
energy  difference  between  the  HOMO  (highest  occupied  molecular  orbital)  and  LUMO  (lowest 
unoccupied  molecular  orbital).  Therefore,  the  study  of  the  electronic  properties  including  the  bandgaps 
may  give  a  clue  of  the  impact  sensitivity  under  hydrostatic  compression.  The  bandgap  of  P-HMX  as  a 
function  of  hydrostatic  pressure  is  illustrated  in  the  Figure  2.7.  It  shows  that  the  bandgap  monotonically 
decreases  with  increasing  hydrostatic  pressure.  This  trend  could  be  understood  as  the  decrease  of 
intermolecular  space  under  compression,  which  leads  to  an  increase  of  different  groups  of  electronic 
bands  and,  hence,  to  the  increase  of  charge  overlap  and  delocalization. 

It  is  worth  mentioning  that  the  rate  of  decrease  of  the  bandgap  is  not  homogeneous:  at  lower  pressure 
up  to  12GPa,  the  rate  is  much  larger  than  those  beyond  12  GPa.  This  indicates  that  the  bandgap  energy 
reduction  is  more  pronounced  in  the  low-pressure  region  than  high-pressure  region. 
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Figure  2.7  The  electronic  energy  bandgap  of  /?-HMX  as  a  function  of  hydrostatic  pressure. 

With  the  monotonic  decrease  of  the  bandgap,  one  could  expect  that  there  is  a  bandgap  closure  where 
HMX  becomes  metallic.  It  is  also  proposed  that  the  compression  in  the  detonation  fronts  causes  local 
metalization  which  attributed  to  the  shear  strains  [95].  The  bandgap  closing  produces  distortion-induced 
molecular  electronic  degeneracy  of  the  highest  occupied  and  lowest  unoccupied  molecular  orbitals  of  an 
energetic  molecule  [96].  The  linear  fitting  of  the  bandgaps  in  range  of  12-100  GPa  gives  the  reduction 
ratio  of  -0.0084  eV/GPa.  The  extrapolation  of  such  a  linear  relationship  implies  that  the  zero  bandgap 
occurs  at  296  GPa.  Our  result  is  consistent  with  a  previous  theoretical  prediction  of  320  GPa  [78]. 

Studies  on  the  exitonic  mechanism  of  detonation  show  that  the  reduction  of  the  bandgaps  promote 
the  HOMO-LUMO  transition  within  a  molecule  [96,  97],  thus  increases  the  sensitivity.  The  smaller  the 
bandgap  is,  the  higher  the  possibility  that  the  electron  transfers  from  valence  band  to  conduction  band, 
causing  decomposition  that  may  lead  to  explosion.  This  correlation  between  the  bandgap  and  detonation 
sensitivity  could  be  understood  as  follows.  The  number  of  excited  sates  increases  as  a  response  to  the 
reduction  of  the  bandgap.  The  increment  of  the  population  of  excited  states  augments  the  possibility  of 
the  chemical  reactions.  As  a  consequence,  it  enhances  the  sensitivity  [98].  Therefore,  the  monotonic 
decrease  of  the  electronic  energy  bandgap  in  Figure  2.7  indicates  that  the  impact  sensitivity  increases  with 
an  increase  of  hydrostatic  pressure. 

2.6  Equation  of  state  of  HMX  at  finite  temperature  using  AIMD  modeling 

Thermal  dynamics  play  a  crucial  role  in  the  behavior  of  energetic  materials.  Especially,  the  thermal 
dynamics  at  high  temperature  could  cause  the  decomposition  of  the  energetic  materials.  Accurate 
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modeling  of  the  energetic  materials  at  finite  temperature  using  first-principles  study  is  computationally 
demanding  and  challenging,  since  it  is  an  extension  of  the  standard  first-principles  calculation  at  zero 
Kelvin  temperature. 

We  studied  the  equation  of  states  of  P-HMX  at  finite  temperature  using  ab  initio  molecular  dynamics 
(AIMD)  with  the  VASP  package.  The  simulation  box  contains  8  HMX  molecules,  which  has  224  atoms. 
The  simulation  box  is  initially  compressed,  to  make  the  P-HMX  at  different  densities.  There  are  four 
densities  we  studied  here,  1.90,  2.18,  2.42,  3.05  g/cmA3,  corresponding  to  the  pressure  of  0,  4,  10,  40  GPa 
at  zero  temperature,  respectively.  At  each  density,  the  system  was  firstly  fully  relaxed,  including  the  cell 
and  ions.  Then  the  system  is  subjected  to  the  AIMD  modeling  with  canonical  ensemble  (NVT).  The  time 
step  is  1.0  femtosecond  (1015  s).  Eight  temperatures  ranged  from  50K  to  400K  are  studied.  These 
temperatures  were  selected  to  study  because  they  are  below  the  decomposition  temperature  of  P-HMX, 
which  is  around  440K.  The  total  time  of  the  simulation  of  each  NVT  ensemble  is  5.0  picoseconds  (10"12s). 


Figure  2.8  Equation  of  state  of  P-HMX  at  different  densities  and  temperatures. 


The  results  of  the  pressure  as  a  function  of  temperature  are  plotted  in  Figure  2.8  for  the  4  densities.  It 
shows  that  the  pressure  linearly  increases  with  temperature.  It  should  be  noted  that  the  canonical 
ensemble  (NVT)  was  used  for  each  state  (point)  using  AIMD  modeling  for  5  picoseconds. 

2.7  Mechanical  response  of  p-HMX  at  the  continuum  level 

We  report  a  theoretical  prediction  of  the  mechanical  properties  of  the  p-HMX  energetic  molecular 
crystal  using  DFT-D3,  a  first-principles  calculation  based  on  density  functional  theory  (DFT)  with  van 
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der  Waals  (vdW)  corrections.  The  mechanical  respons  of  P  -HMX  at  the  continuum  level  corresponding 
to  the  elastic  constants  predicted  from  DFT-D3  calculations  are  compared  with  experiments,  DFT-D2  and 
standard  DFT  calculations,  and  other  theoretical  calculations,  as  in  Figure  2.9. 


Expt  (Rae  et  al.,  2006) 

Expt.  (Sun  etalL,  2009} 

Expt.  {Stevens  and  Eckhardt,  2005) 
^4- Expt.  (Sewell  et  al.,  2002) 

■•Expt.  (Zaug,  1998} 
r-DFT 
-$-DFT-D2 
-8-DFT-D3 

-*-MD  (Sewell  et  al.,  2003} 

-i-MD  (Cui  et  al.,  2010) 

DFT  (Chen  et  al.,  2011) 


Figure  2.9  Mechanical  response  of  P-HMX. 


In  Figure  2.9,  the  mechanical  response  of  P-HMX  at  the  continuum  level  corresponding  to  the  elastic 
constants  predicted  from  DFT-D3  calculations  are  compared  with  experiments,  DFT-D2  and  standard 
DFT  calculations,  and  other  theoretical  calculations.  The  squared  lines  are  the  direct  experimental 
measurements,  opposite  to  others  as  interpolation  from  elastic  constants.  The  uniaxial  compressive  strains 
are  applied  (a)  on  the  (110)  plane,  and  (b)  along  the  [001]  direction,  (c)  The  predicted  behavior  of  uniaxial 
compression  along  the  [001]  direction  from  standard  DFT  calculations  compared  with  the  DFT-D2  and 
DFT-D3  calculations,  and  experiment. 

Our  results  show  that  the  DFT-D3  method  models  the  mechanical  response  of  P-HMX  better  than 
DFT-D2  method.  The  standard  DFT  modeling  without  van  der  Waals  corrections  are  not  acceptable,  as 
demonstrated  in  Figure  2.9(c).  Our  study  suggested  that  the  van  der  Waals  interactions  are  critically 
important  t  in  modeling  the  mechanical  properties  of  this  molecular  crystal. 
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Chapter  3-  QM  study  of  the  mechanical  properties  and  the  equation  of  state  of 

Aluminum  and  Aluminum  oxide 

It  is  desirable  to  further  increase  the  power  and  thrust  output  after  explosion  of  explosives.  One  of  the 
approach  is  the  addition  of  metal  particles  to  energetic  materials.  Among  them,  Aluminum(Al)  powders 
are  widely  used  in  pyrotechnics,  rocket  propellants,  and  explosives.  The  aluminized  explosives  are 
widely  used  energetic  materials  in  the  military  explosives  to  increase  the  air  blast,  incendiary  effects,  and 
bubble  energies  in  underwater  weapons.  A1  is  also  commonly  added  to  propellants  to  increase  the  thrust. 
There  is  an  increased  desire  to  develop  combined-effects  aluminized  explosives  to  achieve  both  excellent 
metal  pushing  and  high  blast  energies  [99],  including  the  study  of  the  detonation  characteristics  of 
aluminized  explosives  [100-103],  to  understand  and  predict  the  macroscopic  response  of  energetic 
materials  on  the  basis  of  mechanical  properties,  which  are  required  in  the  formulation  and 
parameterization  of  mesoscale  continuum  mechanical  constitutive  models  [103-111].  The  accurate 
prediction  of  the  mechanical  properties  and  crystallographic  parameters  are  imperative  to  the 
development  of  physically  based  constitutive  models  that  are  necessary  for  continuum  scale  modeling  of 
deformation  mechanisms,  which  affect  the  sensitivity  of  energetic  molecular  crystals  such  as  HMX. 

The  aluminum  powder  added  into  the  mixture  improves  the  heat  and  temperature  of  the  explosion  for 
the  mixture  and  increases  the  action  duration  of  the  explosion.  Unlike  the  homogenous  ideal  explosive,  it 
is  difficult  to  predict  the  effect  of  aluminized  explosives  because  its  detonation  is  not  ideal.  The  EoS  for 
detonation  products  is  needed  to  predict  explosion  effects.  The  EoS  for  detonation  products  of  ideal 
homogeneous  explosives  can  be  determined  easily  because  its  detonation  mechanism  is  simpler  and  the 
related  research  is  perfect. 

One  approach  to  obtain  the  EoS  of  a  mixture  is  to  combine  the  EoS  of  each  component.  We  take  the 
aluminized  HMX  as  example.  There  are  two  ingredients  A1  (Aluminum)  and  HMX. 

HMX  is  an  important  secondary  explosive  that  exists  in  three  crystallographic  polymorphs  a,  P,  8,  and 
in  hemi-hydrate  y  form.  The  P21/c  monoclinic  P-phase  molecular  crystal  (Figure  2.1)  is 
thermodynamically  the  most  stable  polymorph  of  HMX  at  room  temperature  [112,16].  P-HMX  is  a  high- 
density  energetic  material  with  a  high  detonation  velocity,  which  is  most  commonly  used  in  military  and 
industrial  applications  including  polymer-  bonded  explosives  and  propellant  formulations.  Due  to  its 
importance,  P-HMX  has  been  the  subject  of  extensive  studies  both  experimental  (for  example,  Ref.  [73, 
86,  88])  and  theoretical,  which  includes  constitutive  models  [104-111]  molecular  dynamics  simulations 
[87,  113-115]  and  first-principles  calculations  [45,  49,  82].  The  structures,  mechanical  properties, 
equations  of  state,  and  electronic  properties  of  P-HMX  under  hydro  static  pressure  were  studied  using  a 
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DFT-D2  method  [23].  The  predictions  of  elastic  constants  under  various  functionals  were  reported,  as 
well  as  EOS. 

Aluminum  is  a  common  metal  used  in  everyday  life.  The  advantage  of  using  aluminum  as  an 
ingredient  in  aluminized  explosives  is  because  of  the  high  performance,  as  A1  is  a  very  active  metal  with 
high  energy  density.  Another  advantage  of  aluminized  energetic  materials  is  that  their  reaction  products 
are  environmentally  friendly.  Finally,  the  materials  have  in  general  relatively  low  production  cost.  To 
make  our  study  complete  and  consistent,  we  re-investigate  the  properties  and  EoS  of  Aluminum  using 
newly  developed  more  accurate  Density  Functional  Theory  calculations. 

Mechanical  properties  including  elastic  constants  are  among  the  most  important  physical 
characteristics  of  solid  energetic  materials  as  they  provide  a  description  of  the  mechanical  behavior  at  the 
continuum  level.  Although  first-principles  calculations  based  on  density  functional  theory  (DFT)  provide 
overall  better  predictive  power  than  pure  force-field  models  (for  example,  in  Ref.  [116]),  they  fail  in 
predicting  lattice  constants,  volumes,  and  elastic  constants  of  energetic  materials  with  standard 
approximations,  partially  due  to  their  poor  descriptions  of  dispersion  forces  in  molecular  crystals  [23]. 

We  consider  a  conventional  unit  cell  containing  4  atoms  with  periodic  boundary  conditions  for  face- 
centered  cubic  structure  which  is  the  most  stable  phase  for  Aluminum  under  ambient  condition.  The  total 
energies  of  the  system,  forces  on  each  atom,  stresses,  and  stress-strain  relationships  of  A1  under  the 
desired  deformation  configurations  are  characterized  via  first-principles  calculations  based  on  density- 
functional  theory  (DFT).  DFT  calculations  were  carried  out  with  the  Vienna  Ab-initio  Simulation 
Package  (V ASP) [66,  67]  which  is  based  on  the  Kohn-Sham  Density  Functional  Theory  (KS-DFT) [68-69] 
with  the  generalized  gradient  approximations  as  parameterized  by  Perdew,  Burke,  and  Ernzerhof  (PBE) 
for  exchange-correlation  functions [70]. The  electrons  explicitly  included  in  the  calculations  are  the  3s23pl 
for  aluminum  atoms,  and  2s22p6  for  oxygen  atoms.  The  core  electrons  are  replaced  by  the  projector 
augmented  wave  (PAW)  and  pseudo-potential  approach  [71,  72].  The  kinetic-energy  cutoff  for  the  plane- 
wave  basis  was  selected  to  be  800  eV  in  this  study.  The  calculations  are  performed  at  zero  temperature. 
We  applied  the  PBEsol  functionals  in  this  study. 

The  criterion  to  stop  the  relaxation  of  the  electronic  degrees  of  freedom  is  set  by  total  energy  change  to 
be  smaller  than  0.000001  eV.  The  optimized  atomic  geometry  was  achieved  through  minimizing 
Hellmann-Feynman  forces  acting  on  each  atom  until  the  maximum  forces  on  the  ions  were  smaller  than 
0.001  eV/A.  The  atomic  structures  of  all  the  deformed  and  undeformed  configurations  were  obtained  by 
fully  relaxing  the  unit  cell.  The  simulation  invokes  periodic  boundary  conditions  for  the  three  directions. 
The  irreducible  Brillouin  Zone  was  sampled  with  al6x  16x16  Gamma-centered  k- mesh.  The  initial 
charge  densities  were  taken  as  a  superposition  of  atomic  charge  densities. 
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We  firstly  optimized  the  geometry  of  the  simulation  box  and  obtained  the  lattice  constant  of  4.011  A, 
which  agrees  well  with  experiment  4.05  A  and  previous  DFT  studies  ranged  from  3.97-4.10  A. 

We  examined  the  electronic  band  structure,  which  is  illustrated  in  the  Figure  3.1.  The  Fermi  energy  is 
set  to  be  0.0  eV.  The  T,  M,  X,K  are  the  four  highest  symmetric  points  in  the  first  Brillouin  zone,  as 
(0,0,0), (0.5, 0,0), (0.5, 0.5,0),  and  (0.5, 0.5, 0.5)  respectively.  The  energy  ranged  from  -10.3  eV  to  16.5  eV. 
The  electron  Density  of  State  (DOS)  are  also  plotted  aside  to  the  electronic  band  structure.  It  clearly 
shows  the  consistency  with  each  other. 


Figure  3.1  The  electronic  band  structure  of  FCC  aluminum  (left)  and  the  corresponding  electron 

Density  of  State  (right). 

In  Figure  3.1,  T,  M,  X,  K  are  the  four  highest  symmetric  points  in  the  first  Brillouin  zone,  as 
(0,0,0), (0.5, 0,0), (0.5, 0.5,0),  and  (0.5, 0.5, 0.5)  respectively.  The  Fermi  energy  are  set  to  be  zero. 

To  gain  more  insight  about  the  property-structure  relationships,  we  also  studied  the  partial  electronic 
density  of  state,  as  illustrated  in  Figure  3.2.  It  clearly  shows  that  the  electrons  on  Fermi  surface  are  mainly 
in  p  orbitals.  Further  study  of  the  components  of  the  p  electrons  of  A1  reveals  that  the  pz  electrons  are 
dominant,  which  could  be  the  key  for  its  active  chemical  response. 
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Figure  3.2:  The  projected  electron  density  of  state  of  Aluminum: 

line)  and  the  p  orbital  (red  line). 


Total  (gray  line),  the  s  orbital  (blue 
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Figure  3.3  The  projected  electron  density  of  state  of  p  electrons  in  FCC  Aluminum:  Total  p  electrons 
and  it  its  three  components  along  x,y,z  directions. 
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We  next  studied  the  equation  of  state  of  FCC  aluminum.  We  isotopically  apply  external  pressure  on 
the  simulation  box.  The  corresponding  energy  and  volume  after  fully  relaxation  then  calculated  from  DFT. 
The  internal  energy  as  a  function  of  volume  is  plotted  in  Figure  3.3.  In  addition,  we  fit  the  E-V  curves 
using  the  third-order  Birch-Murnaghan  (BM3)  isothermal  equation  of  state,  namely, 
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Eqn.  3.1 


E(V)  =  E0  + 


9  V0B0 


16 


m2,s  -  > 


3B'0  + 


®2/3-'  *«-«(*) 


2/3 


The  fitting  parameters  were  the  ground  state  free  energy  E0,  the  unit  cell  volume  per  atom  at  zero 
pressure  V0 ,  bulk  modulus  B0 ,  and  £?q — the  pressure  derivative  of  the  bulk  modulus  at  zero  pressure.  Our 
results  are  V0=1632  A3,  bulk  modules  B0=  75  GPa,  B'0=434 ,  and  E0  =  —3.78  eV  per  atom.  Our  results 
have  good  agreement  with  experimental  values. 

The  unit  cell  volumes  calculated  at  different  pressures  were  used  to  generate  P-V  curves,  which  were 
fit  with  the  third-order  Birch-Murnaghan  isothermal  equation  of  state  for  pressure,  namely, 


POO  =  3-y 


7/3  /Kn\5/3 


l  +  f(*o-4) 


[&)-&) 

The  results  of  P-V  curves  and  the  fitting  are  shown  in  Figure  3.4. 


Eqn.  3.2 


Figure  3.4  The  free  energy  as  a  function  of  volume  of  Aluminum  from  first-principles  calculations 
(red  dots)  are  fitted  with  the  third-order  Birch-Murnaghan  (BM3)  isothermal  equation  of  state. 

In  Figure  3.4,  the  fitting  parameters  are  the  ground  state  free  energy  E0,  the  unit  cell  volume  per  atom 
at  zero  pressure  V0,  bulk  modulus  B0 ,  and  B'0 — the  pressure  derivative  of  the  bulk  modulus  at  zero 
pressure. 

In  Figure  3.5,  the  fitting  parameters  were  the  unit  cell  volume  per  atom  at  zero  pressure  V0 ,  bulk 
modulus  B0 ,  and  B'0 — the  pressure  derivative  of  the  bulk  modulus  at  zero  pressure.  The  dashed  vertical 
line  denotes  the  position  of  unstrained  state. 
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Under  ambient  condition,  the  A1  powder  or  nanoparticles  will  have  chemical  reactions  with  oxygen  in 
air  and  generate  a  thin  layer  of  A1203  on  the  surface.  This  is  a  protective  layer  which  prevent  further 
oxidation  of  Al.  So  A1203  could  be  a  considerable  ingredient  when  the  size  of  the  A1  powder  is  in 
nanometers.  The  knowledge  of  the  mechanical  properties,  thermal  properties,  and  chemical  reactions 
related  to  A1203  are  critical  important  for  the  understanding,  design,  and  engineering  the  aluminized 
explosives.  To  that  end,  we  have  studied  the  structure,  elastic  properties,  and  electronic  properties. 


Figure  3.5:  The  pressure  as  a  function  of  volumetric  strain  of  aluminum  from  first-principles 
calculations  (red  dots)  are  fitted  with  the  third-order  Birch-Murnaghan  (BM3)  isothermal  equation  of  state. 

The  alpha  A1203  has  a  crystal  structure  of  R-3c  with  space  group  number  of  167.  The  conventional 
unit  cell  contains  12  Al  atoms  and  18  oxygen  atoms.  The  full  optimization  of  the  geometry  using  PBEsol 
functionals  reveals  the  lattice  constants  of  a  =  4.777  A  and  c  =  13.024  A,  agree  well  with  experiment  of 
a  =  4.758  A  and  c  =  12.991  A,  respectively.  The  simulation  cell  with  the  atoms  is  displayed  in  Figure 
3.6.  The  grey  atoms  are  Al  and  the  red  atoms  are  oxygen. 

In  Figure  3.6,  the  crystal  structure  is  R-3c  with  space  group  number  of  167.  The  lattice  constants  are 
a  =  4.777  A  and  c  =  13.024  A  predicted  from  PBEsol  functions. 

Then  we  isotopically  apply  external  pressure  on  the  simulation  box.  The  corresponding  energy  and 
volume  after  fully  relaxation  then  calculated  from  DFT.  The  internal  energy  as  a  function  of  volume  is 
fitted  the  third-order  Birch-Murnaghan  (BM3)  isothermal  equation  of  state  (Eqn.  3.1).  The  results  are 
illustrated  in  Figure  3.7. 
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In  Figure  3.7,  the  fitting  parameters  are  the  ground  state  free  energy  E0,  the  unit  cell  volume  per  atom 
at  zero  pressure  V0 ,  bulk  modulus  B0 ,  and  Br0 — the  pressure  derivative  of  the  bulk  modulus  at  zero 
pressure. 


Figure  3.6:  The  conventional  unit  cell  of  the  alpha- A1203  with  12  A1  (grey)  and  18  oxygen  (red)  atoms. 

Our  results  are  V0=8.62  A3  per  atom  (which  is  257.37  8.62  A3  per  unit  cell),  bulk  modules  B0=  224.09 
GPa,  B o=4.19,  and  E0  =  —7.70  eV  per  atom.  The  free  energy  is  -231.0  eV  per  unit  cell.  Our  results  have 
good  agreement  with  experimental  values. 

The  heat  of  formation  of  A1203  can  be  calculated  through  the  energy  difference  between  the  reactants 
and  products,  as  H=E(A1203)-2E(A1)-1.5E(02).  At  T=0K,  the  energy  of  A1  E(Al)=-3.78  eV,  the  energy 
of  A1203  is  -38.5  eV,  the  02  is  -8.855  eV.  The  formation  energy  at  zero  K  is  -17.6575.  Therefore,  the 
heat  of  formation  is  1703.7  kJ/mol. 

The  unit  cell  volumes  calculated  at  different  pressures  were  used  to  generate  P-V  curves,  which  were 
fit  with  the  third-order  Birch-Murnaghan  isothermal  equation  of  state  for  pressure  (Eqn  3.2).  The  results 
are  shown  in  Figure  3.8. 
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Figure  3.7:  The  free  energy  as  a  function  of  volume  of  alpha- A1203  from  first-principles  calculations 
(red  dots)  are  fitted  with  the  third-order  Birch-Murnaghan  (BM3)  isothermal  equation  of  state. 

In  Figure  3.8,  the  fitting  parameters  are  the  unit  cell  volume  per  atom  at  zero  pressure  F0,  bulk 
modulus  BO,  and  BO’  at  zero  pressure.  The  dashed  vertical  line  denotes  the  position  of  unstrained  state. 


Figure  3.8:  The  pressure  as  a  function  of  volumetric  strain  of  alpha-A1203  from  first-principles 
calculations  (red  dots)  are  fitted  with  the  third-order  Birch-Murnaghan  (BM3)  isothermal  equation  of  state 

for  pressure. 
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Chapter  4-Validation  of  the  newly  developed  ReaxFF  potential  ReaxFF-d3 

As  a  molecular  crystal,  the  /?- HMX  need  accurate  energy  formula  to  prescribe  the  van  der  Walls 
interactions  between  molecules  which  is  critical  in  modeling  the  behaviors  of  this  materials.  There  are 
extensive  studies  to  improve  the  modeling  of  dispersion,  or  van  der  Waals  (vdW)  interactions  [68],  such 
as  DFT-D  [51,  55,  57,  60,  61],  vdW-TS  [59,  62],  and  vdW-DF  [52,  63]  methods.  The  DFT-D1  and  DFT- 
D2  methods  are  classified  as  the  1st  rung  of  Jacob’s  ladder  for  DFT-based  dispersion  correction  schemes 
[68].  The  DFT-D3  and  vdW-TS  methods  belong  to  the  second  rung  that  utilizes  environment  dependent 
C6  corrections.  Non-local  density  functionals  like  DFT-DF  and  DFT-DF2  methods  sit  in  the  3rd  rung 
[68],  of  which  overall  accuracy  is  expected  to  increase,  as  well  as  computational  demand.  Recently,  a  few 
improvements  in  vdW-DF  methods  have  been  proposed,  including  optPBE-vdW,  optB88-vdW,  optB86b- 
vdW  [53,  117],  vdW-DF2-C09  [54],  vdW-DF-cx  [118],  and  rev-vdW-DF2  [119].  These  methods  are 
designed  to  minimize  the  error  of  the  binding  energies,  equilibrium  separations,  and  interaction  energy 
curves  in  the  molecular  duplexes  in  a  certain  collection  of  materials,  such  as  the  S22  data  set  [120], 
aiming  for  a  general  use  for  all  materials.  The  accuracy  of  these  models  on  energetic  materials  like  HMX 
have  been  assessed.  With  the  compromise  the  accuracy  and  the  computational  cost,  we  choose  the  vdW- 
D3  as  the  main  functionals  for  out  modeling  of  the  HMX  in  density  functional  theory.  The  mechanical 
properties  and  equations  of  state  were  reported  previously  elsewhere  by  the  authors. 

As  a  second  stage  for  bottom-up  multiscale  modeling,  it  is  vital  to  pass  the  information  from  quantum 
mechanics  electron  level  to  molecular  dynamics  modeling  in  atomic  level  in  a  form  of  force  field.  In  order 
to  capture  the  fast  chemical  reactions  and  kinetics  during  the  detonation,  the  reactive  force  fields  are 
desirable.  To  that  end,  we  adopted  the  ReaxFF  potential  as  our  target.  A  new  set  of  ReaxFF  potential  was 
developed  last  year  with  the  new  training  sets  of  HMX  structures  that  are  produced  from  vdW-D3 
calculations.  As  a  result,  we  denote  this  new  set  of  ReaxFF  potential  for  C-H-N-O  systems  as  ReaxFF- 
D3. 

We  continued  to  validate  ReaxFF-D3  by  its  performance  in  molecular  dynamics  simulations.  We 
studied  the  thermal  expansion  coefficients  of  beta-HMX,  as  the  experimental  reference  is  available. 

We  study  the  volume  change  under  conditions  of  NPT  at  temperature  T  ranged  from  IK  to  500K  and 
the  pressure  of  1  atm  (105  Pa).  The  time  step  in  the  MD  simulations  is  0.25fs.  The  simulations  run  for 
20000  steps  which  is  corresponding  to  5ps.  The  system  contains  14336  atoms  which  is  512  molecules,  or 
256  primitive  unit  cells.  Such  a  large  system  is  used  for  better  statistics  during  MD  simulations. 

The  results  of  our  NPT  simulations  at  varies  temperatures  are  illustrated  in  Figure  4.1.  The  volume  is 
for  primitive  unit  cell  that  contains  56  atoms,  or  2  molecules.  We  observed  that  the  unit  cell  volume 
increases  when  the  temperature  increases.  We  fit  the  results  of  volumes  as  a  function  of  temperature  using 
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a  quadratic  function,  V  =  0.000019734  T2  +  0.074236 T  +  503.254  as  displayed  in  the  Figure  4.1.  The 
thermal  expansion  then  can  be  obtained  by  the  definition. 

dln(y) 


a  = 


dT 


Eqn.  4.1 


The  thermal  expansion  then  is  plotted  in  the  same  figure  to  the  right  axis  (Figure  4.1).  Our  results 
shows  that  the  volume  is  503.23  A3  and  527.28  A3  at  T=0K  and  300K  respectively.  The  thermal 
expansion  is  163.25xl06  at  300K. 


Figure  4.1:  The  volume  as  a  function  of  temperature  of  /?-HMX  modeled  in  molecular  dynamics  using 

ReaxFF-D3  as  potential. 


It  should  be  noted  that  in  Figure  4.1,  the  function  is  fitted  into  a  quadratic  function.  The  thermal 
expansion  is  plotted  to  the  right  axis.  The  original  data  of  volume  is  denoted  by  red  dots. 

As  a  comparison,  we  studied  the  same  systems  but  with  the  ReaxFF-lg  [19],  the  latest  ReaxFF 
potentials  come  with  LAMMPS.  The  NPT  results  of  volumes  as  a  function  of  temperature  are  plotted  in 
Figure  4.2,  including  the  thermal  expansions.  The  ReaxFF.lg  results  show  that  the  volume  is  503.60  A3 
and  527.52  A3  at  T=0K  and  300K  respectively.  The  thermal  expansion  is  165.14xl06  at  300K. 

It  should  be  noted  that  in  Figure  4.2,  the  function  is  fitted  into  a  quadratic  function.  The  thermal 
expansion  is  plotted  to  the  right  axis.  The  original  data  of  volume  is  denoted  by  red  dots. 
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Figure  4.2  The  volume  as  a  function  of  temperature  of  /?- HMX  modeled  in  molecular  dynamics  using 

ReaxFF-lg  as  potential. 


Experimentally,  the  thermal  expansion  of  the  beta-HMX  was  studied  [50]  in  the  temperature  range  of 
-150  to  30  °C.  (3-HMX  is  monoclinic  (a  =  6.5255(10)  A,  b  =  11.0369(18)  A,  c  =  7.3640(12)  A,  and  p  = 
102.67(1)°),  space  group  P21/n.  The  results  of  the  volumes  at  various  temperature  then  plotted  in  the 
Figure  4.3.  The  volumes  are  also  fitted  into  a  quadratic  function.  The  thermal  expansion  then  obtained.  In 
this  way,  we  have  the  experimental  values  of  volume  which  is  500.318  A3  and  516.984  A3  at  T=0  and 
300K,  respectively.  The  thermal  expansion  at  300  is  139.066xl06  at  300K. 
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Figure  4.3:  The  volume  as  a  function  of  temperature  of  (3- HMX  measured  in  experiment  [50]. 
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It  should  be  noted  that  in  Figure  4.3,  the  volume  is  fitted  into  a  quadratic  function.  The  thermal 
expansion  were  plotted  to  the  right  axis.  The  original  data  of  volume  is  denoted  by  red  dots. 

Besides  the  thermal  expansion,  we  studied  the  vibration  frequencies  of  the  beta-HMX  crystals  at  finite 
temperature.  In  general,  Vibrational  spectroscopy  is  a  sensitive  probe  of  the  atomic  structure  and  of  the 
chemical  bonding  and  thus  of  the  electronic  structure.  In  this  study,  the  phonon  density  of  states  (PDOS) 
has  been  studied  at  T=300K  using  both  ReaxFF-D3  and  ReaxFF-lg  potentials  for  better  comparison.  The 
PDOS  is  calculated  from  the  Fourier  transform  of  the  velocity  autocorrelation  function  [121],  as 


PDOS(ft>)=  ^  /  e-ilBt<£f=1  vjWvjt  0))dt 


Eqn.  4.2 


where  Vj(0)  is  the  average  velocity  vector  of  a  particle  j  at  initial  time,  Vj(t)  is  its  velocity  at  time  t ,  and  co 
is  the  vibration  wave  number.  N  is  the  number  of  atoms  in  the  system.  Our  results  are  plotted  in  the 
Figure  4.4. 


Figure  4.4  The  phonon  density  of  states  of  /?- HMX  at  300K  predicted  by  ReaxFF-D3  (top)  and 

ReaxFF-lg  (bottom)  potentials. 


It  can  be  seen  that  the  main  frequencies  are  represented  and  similar  in  both  potentials.  For  example, 
the  largest  frequency  goes  up  to  97.2  THz  or  3243  cm"1.  This  is  the  typical  vibration  of  hybridized  C-H 
bonds.  The  nitro  groups  (-N02)  has  vibration  frequency  of  1550  cm"1  or  46.5  THz  are  also  captured. 
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Chapter  5-QM  study  of  the  stabilities,  mechanical  properties  and  equations  of 

state  of  Iron  under  extreme  conditions 


Iron  is  one  of  the  most  important  metals  in  everyday  life.  It  is  the  most  common  element  in  the  whole 
planet  Earth,  forming  much  of  Earth’s  outer  and  inner  core.  Iron  is  the  mostly  used  metal  case  for 
explosives.  To  have  better  description  of  the  performance  of  detonation  within  a  metal  case,  we  need  to 
have  better  understand  of  the  stabilities,  mechanical  properties,  thermal  properties,  and  equations  of  state 
of  Iron  under  extreme  conditions. 


The  thermal  properties  of  solids  at  constant  volume  can  be  calculated  from  their  phonon  density  of 
states  as  a  func-  tion  of  frequencies.  The  phonon  contribution  to  the  Helmholtz  free  energy  F  A  is  given 
by: 

hteq/v 

FA  =  2'Lq,vhUq,v  +  kBTYJq,v\n[l-e  kBT]  Eqn.  5.1 

where  kB  and  h  are  the  Boltzmann  constant  and  the  reduced  Planck  constant,  respectively,  q  and  v  are  the 
wave  vector  and  band  index,  respectively,  ooq9V  is  the  phonon  frequency  at  q  and  v,  and  T  is  the 
temperature. 

The  heat  capacity  at  constant  volume  Cv  and  the  entropy  S  at  constant  volume  are  given  by  the 
following  formula: 


Cy=^q,vkBCp2  + 

ha)qiV 

e  kBT 

Eqn.  5.2 

hcoqv 

(e  kBT  -1)2 

S  =  -kBY.q,v ln  1  -  e 

hu>q)V~\ 

kBT  ___ 

h^q,v 

Eqn.  5.3 
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e 
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In  practical,  thermodynamical  problems  related  to  solids,  the  thermal  properties  need  to  be  known  at 
constant  pressure.  They  can  be  calculated  from  the  previous  quantities  through  thermodynamic 
relationship.  The  Gibbs  free  energy  G  may  be  written  as 

G(T,P)  =  min v[U(V)  +  FA(F,V )  +  PV]  Eqn.  5.4 

where  V  and  P  are  the  volume  and  the  pressure,  respectively.  U(V)  is  the  total  energy  of  electronic 
structure  at  constant  volume.  The  right-hand  side  of  Eq.  4  means  the  function  inside  the  square  brackets  is 
minimized  with  respect  to  the  volume  at  each  couple  of  T  and  p  variables. 

The  heat  capacity  at  constant  pressure  Cp  is  derived  from  G(T,  p)  by 


d2G(T,P ) 


,  dV(T,P)  dS(T,V) 


\v=v(t,p)  +  CV(T,V(T,Py)  Eqn.  5.5 


where  V(T,  P)  is  the  equilibrium  volume  at  temperature  T  and  P. 
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We  have  employed  quasiharmonic  approximation  (QHA)  to  calculate  the  thermal  properties  at 
constant  pressure.  U(V)  and  FA(T,V )  were  calculated  at  11  volume  points,  from  -15%  to  +15%  of 
volume  change  around  the  equilibrium.  The  corresponding  pressure  is  ranged  from  42.6  GPa  to  -17.9  GPa. 

The  thermodynamic  functions  were  fitted  to  the  integral  form  of  the  Vinet  equation  of  state  [99]  (EOS) 
at  p  =  0.  Gibbs  free  energies  at  finite  temperatures  were  obtained  as  the  minimum  values  of  the 
thermodynamic  functions,  and  the  corresponding  equilibrium  volumes  and  isothermal  bulk  moduli  were 
obtained  simultaneously  from  the  Vinet  EOS.  Unit  cells  used  to  calculate  U(V)  and  FA(T,V )  were 
relaxed  by  the  first-principles  calculation  under  the  hydrostatic-stress  conditions.  These  procedures 
applied  for  bee  iron  are  demonstrated  in  Figure  5.1,  where  U(V)  +  FA(T,V ))  as  a  function  of  unit-cell 
volume  at  temperatures  are  shown. 

The  CP  of  bee  iron  as  a  function  of  temperature  obtained  from  the  quasi-harmonic  approximation  are 
shown  in  Figure  5.2,  ranging  from  0  to  2000K.  Please  keep  in  mind  that  the  melting  temperature  of  iron  is 
1800K.  The  crystal  structure  might  not  be  stable. 

The  free  energy,  entropy,  and  heat  capacity  at  constant  volume  Cv  as  a  function  of  temperature  of 
BCC  iron  from  QHA  are  obtained  and  plotted  in  Figure  5.3. 

The  phonon  band  structure  of  BCC  iron  under  pressure  of  0  GPa  and  43  GPa  are  plotted  in  Figure  5.4. 
There  is  negative  frequency  at  high  pressure  of  43  GPa,  which  indicates  the  instability  of  BCC  iron  under 
such  pressure.  In  fact,  there  is  a  phase  change  from  BCC  structure  to  HCP  structure.  At  pressures  above 
approximately  10  GPa  and  temperatures  of  a  few  hundred  kelvin  or  less,  a-iron  changes  into  a  hexagonal 
close-packed  (hep)  structure,  which  is  also  known  as  s-iron  or  hexaferrum. 


Figure  5.1:  U(V)  +  FA(T,  V)  as  a  function  of  unit-cell  volume  of  bee  iron. 
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In  Figure  5.1,  the  solid  curves  show  the  fitted  thermodynamic  functions.  The  minimum  values  of  the 
fitted  thermodynamic  functions  at  temperatures  are  depicted  by  the  crosses.  The  dashed  curve  passing 
through  the  crosses  is  guide  to  the  eye. 


Figure  5.2  The  CP  of  bee  iron  as  a  function  of  temperature  obtained  from  the  quasi-harmonic 


approximation. 


Figure  5.3  The  free  energy,  entropy,  and  heat  capacity  at  constant  volume  Cv  as  a  function  of 

temperature  of  BCC  iron  from  QHA. 
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Figure  5.4  Phonon  band  structure  of  BCC  iron  under  pressure  of  0  GPa  (top)  and  43  GPa  (bottom). 


There  is  negative  frequency  at  high  pressure  of  43  GPa,  which  indicates  the  instability  of  BCC  iron 
under  such  pressure. 
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Chapter  6-Molecular  dynamics  (MD)  simulation  using  LAMMPS 

6.1  Hugoniot  EoS  of  solid  HMX 

The  Hugoniot  EoS  describes  the  locus  of  thermodynamic  states  of  solid  high  explosives  after  a  shock 
passing  through.  It  can  be  derived  from  the  conservation  of  mass,  momentum,  and  energy,  known  as 
Rankine-Hugoniot  conditions  [122]: 

Pius  —  P2(^5  —  ut)  Conservation  of  mass  Eqn.  6.1 

P2  ~Pi  ~  P2U2(VS  —  u2 )  =  p\Usu2  Conservation  of  momentum  Eqn.  6.2 
p2u  2  —  Pius(.2Uz  +  E2  —  E  i)  Conservation  of  energy  Eqn.  6.3 

where  us  is  the  shock  wave  speed,  p2  and  p1  are  the  density  of  material  behind  and  in  front  of  the  shock, 
u2  is  the  particle  velocity  of  the  material  behind  the  shock,  p2  and  are  the  pressure  of  the  material 
behind  and  in  front  of  the  shock,  and  E2  and  E1  are  the  specific  internal  energies  in  the  two  regions. 

Based  on  Rankine-Hugoniot  conditions,  the  following  relations  can  be  obtained: 


U2_. 

Eqn.  6.4 

US  V± 

_  \V2~Vl_ 

—  ^1 

Eqn.  6.5 

vl~v2 

U2  =  V(P2  -Pl)Oi  -  V2) 

Eqn.  6.6 

El=\  (P2  +  Pi)  (“““)=  I  (P2  +  Pl)Ol  -  V2) 

Eqn.  6.7 

where  v2  and  v1  are  the  specific  volumes  of  the  material  behind  and  in  front  of  the  shock. 

The  parameters  of  the  JWL  EoS  of  solid  PBX  9501  are  listed  in  Table  6.1  [24].  PBX  9501  is  chosen 
because  it  consists  of  95  weight  %  HMX,  2.5  weight  %  estane  binder,  and  2.5  weight  %  BDNPA/F. 


Table  6.1  Parameters  of  JWL  EoS  of  solid  PBX  9501  [24] 


A  ( GPa ) 

B  (GPa) 

ri 

r2 

o>  Cv(J/(Kg  ■  K)) 

T0(K ) 

732000 

-5.2654 

14.1 

1.41 

0.8867  2.7806  x  106 

298 

6.2  JWL  EoS  of  gaseous  products  of  HMX 

The  expansion  of  detonation  products  of  high  explosives  is  considered  an  isentropic  process.  It 
includes  complicated  chemical  reactions,  occurring  in  very  short  time.  Cylinder  test  has  long  been  the 
principal  method  to  obtain  the  data  needed  for  fitting  parameters  of  the  JWL  EoS  of  explosive  detonation 
products  [17].  It  consists  of  detonating  a  cylinder  of  explosive  confined  by  copper  and  measuring  the 
velocity  of  the  expanding  copper  wall.  The  standard  tube  has  a  1-inch  internal  diameter,  0.1  in  copper 
wall,  and  is  12  in  long.  The  configuration  of  cylinder  test  is  shown  in  Figure  6.1. 
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Detonator 


6.1  Configuration  of  cylinder  test 

The  parameters  of  the  JWL  EoS  of  gaseous  products  of  PBX  9501  are  listed  in  Table  6.2  [24]. 
Table  6.2  Parameters  of  JWL  EoS  of  gaseous  products  of  PBX  9501 


A  ( GPa ) 

B  (GPa) 

n 

r2 

CO 

Cv(J/(Kg-K )) 

T0(K) 

1668.9 

59.69 

5.9 

2.1 

0.45 

1.0  x  106 

298 

6.3  Structure  and  relaxation  of  p-HMX  using  LAMMPS 

We  focus  on  the  P-HMX  (octahydro-l,3,5,7-tetranitro-l,3,5,7-tetrazocine),  a  typical  secondary 
explosive  with  chemical  formula  of  C4H8N808.  The  geometry  of  the  primitive  unit  cell  of  the  P-HMX  is 
shown  in  Figure  6.2(a)  with  the  cell  box,  which  is  obtained  from  our  previous  DFT-D2  study  [23].  The  P- 
HMX  is  a  P21/c  monoclinic  crystal  with  space  group  number  14.  The  lattice  parameters  are  listed  in 
Table  2.1.  Replication  of  the  primitive  unit  cell  of  HMX  is  shown  in  Figure  6.2(b),  (a),  where  C,  H,  O, 
and  N  atoms  are  represented  by  yellow,  aqua,  red,  and  grey  balls,  respectively.  We  tested  three  different 
systems  with  4x2x4,  5x3x5,  and  8x5x7  primitive  unit  cells,  with  1792,  4200,  and  15680  atoms, 
respectively  and  it  is  found  that  the  results  converge  at  the  system  with  1792  atoms.  Therefore,  we  report 
our  results  from  the  simulation  cells  with  1792  atoms  as  depicted  in  Figure  6.2(b). 
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Figure  6.2  Simulation  cells  and  geometry  of  [3-HMX 
The  structure  of  |3-HMX  at  OK  was  obtained  from  our  previous  DFT-D2  calculations  [23].  The  system 
which  contains  1792  atoms  is  relaxed  using  the  NPT  (constant  atom  number  N,  pressure  P ,  temperature 
T ),  NVT  (constant  atom  number  N,  volume  V,  temperature  T),  and  NVE  (constant  atom  number  N, 
volume  V,  energy  E )  ensembles  successively  with  ReaxFF-d3.  Due  to  the  inadequate  consideration  of  van 
der  Waals  (vdW)  attractions,  when  using  original  ReaxFF  to  relax  the  solid  molecule,  the  equilibrium 
volume  will  be  10-15%  higher  than  experiment  [19].  This  disadvantage  is  overcome  by  ReaxFF-lg  (and 
ReaxFF-d3),  which  gives  sufficient  consideration  of  vdW  interactions  in  its  description  of  reactive  force 
fields.  In  the  structural  of  relaxation,  the  time  step  is  set  to  be  0.25  fs  while  the  simulation  time  is  set  to 
be  50  ps  for  each  ensemble. 

The  lattice  constants  for  the  unit  cell  of  [3-HMX  from  experiment  [88],  DFT-D2  [23]  studies  are  listed 
in  Table  6.3. 


Table  6.3  Lattice  constants  of  [3-HMX  from  experiment  [88],  DFT-D2  calculations  [23] 


a  (A) 

b(  A) 

c(A) 

a 

P 

Y 

Experiment[88] 

6.54 

11.05 

8.70 

90 

124.30 

90 

DFT-D2[23] 

6.542 

10.842 

8.745 

90 

124.41 

90 

The  supercell  of  [3-HMX  consists  of  64  |3-HMX  molecules,  which  contains  1792  atoms  in  total.  The 
density  of  |3-HMX  before  relaxation  is  1891.8  /m3  .  Three  ensembles  (NPT,  NVT,  NVE)  are  used 
sequentially  to  relax  the  system.  Each  ensemble  takes  50  ps.  The  variation  of  volume,  total  energy, 
temperature  and  pressure  of  the  system  during  relaxation  is  shown  in  Figure  6.3. 
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Figure  6.3  Variation  of  volume  (a),  total  energy  (b),  temperature  (c),  and  pressure  (d)  of  P-HMX 
system  during  the  geometry  relaxation  at  the  temperature  of  300  K  in  MD  simulations. 

It  is  well  known  that  P-HMX  is  very  sensitive  to  fluctuations  in  temperature.  It  will  transform  from  P- 
HMX  to  5-HMX  above  436 K  [123].  Therefore,  it  is  necessary  to  control  temperature  fluctuations  to 
prevent  phase  transformation.  It  can  be  seen  from  Figure  6.3(c)  that,  during  the  relaxation,  the 
temperature  is  limited  to  around  300X,  which  ensures  no  phase  transformation  during  relaxation. 
According  to  Figure  6.3(a),  the  volume  of  the  supercell  fluctuates  significantly  in  the  NPT  stage,  and  the 
maximum  volume  of  the  supercell  reaches  17558.66  A3.  In  the  successive  NVT  and  NVE  stages,  the 
volume  stabilizes  at  16924.73  A3  (or  528.9  A3  per  primitive  unit  cell).  After  relaxation,  the  density  of  P- 
HMX  is  1858.7  kg/m3 ,  which  is  very  close  to  the  experimental  value  (1838.0  kg/m3)  at  ambient 
temperature.  According  to  Figure  6.3(d),  after  the  initial  fluctuation,  the  pressure  of  the  system  fluctuates 
around  0  GPa. 

6.4  Calculation  of  JWL  EoS  of  solid  p-HMX  using  LAMMPS 

The  JWL  EoS  of  solid  P-HMX  is  obtained  through  a  series  of  NVT  simulations  using  LAMMPS  with 
ReaxFF-d3.  NVT  is  preferred  because  as  mentioned  in  section  above,  in  ignition  and  growth  model, 
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temperature  is  introduced  in  the  JWL  EoS  of  solid  explosive.  According  to  Yoo  [16],  the  pressure-volume 
relation  calculated  under  hydrostatic  conditions  is  in  good  agreement  with  Hugoniot  curve.  The 
simulation  scheme  is  as  follows:  compress  the  p-HMX  along  a-direction  of  the  supercell  to  predefined 
volume  and  temperature,  thus,  for  each  set  of  volume  and  temperature,  there  is  a  corresponding  pressure. 
In  the  implementation  of  NVT  simulation,  the  time  step  is  set  to  be  0.1  fs  and  the  total  simulation  time  is 
set  to  be  15  ps.  The  variation  of  volume,  total  energy,  temperature  and  pressure  of  the  supercell  is  shown 
in  Figure  6.4,  where  the  compression  direction  is  a-direction  of  the  supercell  and  the  predefined  volume 
ratio  is  0.5  (compress  the  volume  to  0.5  its  original  volume  V0)  while  the  predefined  temperature  is  300  K. 


Figure  6.4  Variation  of  volume  (a),  total  energy  (b),  temperature  (c),  and  pressure  (d)  of  P-HMX  under 

compression  at  7=1000X,  —  =  0.5. 

^0 

The  comparison  of  the  EoS  calculated  by  MD  with  that  obtained  from  experiment  [16],  and  JWF  EoS 
of  solid  PBX  9501  [24]  is  shown  in  Figure  6.5,  where  x-axis  represents  volume  ratio,  y-axis  represents 
the  pressure  of  the  supercell  of  P-HMX  mentioned  above.  It  should  be  noted  that  the  experiment  was 
performed  at  room  temperature  and  the  JWL  EoS  of  solid  p-HMX  obtained  by  MD  is  calculated  at  300  K. 
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-o-JWL  EoS  of  solid  PBX  9501 


Figure  6.5  Comparison  of  the  EoS  calculated  by  MD  with  that  obtained  from  Yoo’s  experiment  and 

JWL  EoS  of  solid  PBX  9501 

It  can  be  seen  from  Figure  6.5  that  the  EoSs,  either  EoS  obtained  from  MD  simulation  or  EoS  of  solid 
PBX  9501,  do  not  agree  well  with  Yoo’s  experiment.  The  EoS  calculated  by  MD  is  close  to  Yoo’s 
experiment  at  high  volume  ratio  (greater  than  0.97).  However,  it  gradually  deviates  from  experiment 
when  the  system  is  further  compressed.  The  EoS  calculated  by  MD  is  greater  than  EoS  of  solid  PBX  9501 
when  volume  ratio  is  greater  than  0.725.  When  volume  ratio  is  smaller  than  0.725,  the  case  is  opposite.  It 
is  probably  because  that  PBX  9501  contains  95%  HMX  and  5%  other  materials  (whose  density  is  around 
1050  Kg/m3,  smaller  than  the  density  of  HMX,  which  is  1910  Kg/m3) . 

The  parameters  of  JWL  EoS  of  solid  PBX  9501  [24]  and  fitted  parameters  of  EoS  of  solid  P-HMX 
obtained  from  MD  simulation  are  listed  in  Table  6.4. 

Table  6.4  Parameters  of  JWL  EoS  of  solid  PBX  9501  [24]  and  fitted  parameters  of  EoS  of  solid  (3- 

HMX  obtained  from  MD 


As  ( GPa ) 

Bs  (GPa) 

n  s 

r2s 

U)s 

T0(K) 

JWL  EoS 

of  solid  PBX 

732000 

-5.2654 

14.1 

1.41 

0.8867 

298 

9501  [24] 

Fitted 

EoS  of  solid 

129230 

-0.3678 

11.8767 

-2.4892 

4.2431 

300 

p-HMX 
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6.5  Calculation  of  JWL  EoS  of  gaseous  products  of  p-HMX  using  LAMMPS 

Similar  to  the  calculation  of  JWL  EoS  of  solid  p-HMX,  the  JWL  EoS  of  gaseous  products  of  P-HMX 
was  also  obtained  through  NVT  simulations  using  LAMMPS  with  ReaxFF-d3.  Prior  to  the  NVT 
simulation,  NPHUG  is  used  to  compress  the  system  to  around  64  GPa.  NPHUG  performs  a  time 
integration  of  the  Hugoniotstat  equations  of  motion  developed  by  Ravelo  et  al.  [122].  Then,  the  system  is 
dilated  from  compressed  state.  The  simulation  scheme  is:  expand  the  p-HMX  along  a-direction  of  the 
supercell  to  a  predefined  volume  and  temperature,  thus,  for  each  set  of  volume  and  temperature,  there  is  a 
corresponding  pressure.  The  time  step  is  set  to  be  0.1  fs  and  the  total  simulation  time  is  set  to  be  15  ps. 
The  variation  of  volume,  total  energy,  temperature  and  pressure  of  the  supercell  is  shown  in  Figure  6.6, 
where  the  predefined  volume  ratio  is  4  (expand  the  volume  to  4  times  its  compressed  volume  Vc ),  and  the 
predefined  temperature  is  3000X. 
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Figure  6.6  Variation  of  volume  (a),  total  energy  (b),  temperature  (c),  and  pressure  (d)  of  P-HMX  in 

NVT  expansion,  T=1000X,  ^  =  4.0. 
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The  comparison  of  JWL  EoS  of  gaseous  products  of  P-HMX  obtained  through  MD  with  JWL  EoS  of 
PBX  9501  is  shown  in  Figure  6.7,  where  x-axis  represents  volume  ratio,  y-axis  represents  the  pressure  of 
the  supercell  of  P-HMX  mentioned  above.  Temperature  is  set  to  be  3000  X,  which  is  close  to  the 
temperature  of  PBX  9501  at  CJ  state. 


Figure  6.7  Comparison  of  JWL  EoS  of  gaseous  products  of  P-HMX  with  JWL  EoS  of  PBX  9501. 

It  can  be  seen  from  Figure  6.7  that  the  EoS  calculated  via  MD  agrees  very  well  with  JWL  EoS  of 
gaseous  products  of  PBX  9501  in  ignition  and  growth  model. 

The  parameters  of  the  JWL  EoS  of  gaseous  products  of  PBX  9501  [24]  and  fitted  parameters  of  EoS 
of  /?-HMX  obtained  from  MD  simulation  are  listed  in  Table  6.5. 

Table  6.5  Parameters  of  the  JWL  EoS  of  gaseous  products  of  PBX  9501  [24]  and  fitted  parameters  of 


EoS  of  jff-HMX  obtained  from  molecular  dynamics 


Ag  (GP  cl) 

Bg  (GP  cl) 

rl  9 

r2  9 

COg 

To  VO 

EoS  of  gaseous 

products  of  PBX 

1668.9 

59.69 

5.9 

2.1 

0.45 

3000 

9501 

Fitted 


parameters  of  EoS 

1015  85.1  5.35  2.53  0.7041  3000 

of  gaseous  products 
of  jff-HMX 
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Chapter  7-QM-MD-SPH  simulation  of  energetic  materials 

7.1  Shock  wave  and  detonation  model 

Generally,  explosive  can  be  categorized  into  two  classes  [124]:  low  and  high  explosive.  Low 
explosives  are  often  used  as  propellants  in  rocket  motors  and  fireworks.  When  detonated,  the  chemical 
reaction  in  low  explosive  will  produce  a  sub-sonic  deflagration  wave.  In  contrast  to  low  explosive,  the 
chemical  reaction  in  high  explosive  will  generate  a  super-sonic  detonation  wave.  According  to  the 
constituents,  high  explosive  can  be  further  classified  to  homogeneous  explosive  and  heterogeneous 
explosive.  Homogeneous  explosive  has  a  uniform  composition,  such  as  liquid  nitromethane  [125]. 
Heterogeneous  explosive  consists  of  poly  crystalline  materials  containing  voids  of  various  shapes  and 
sizes,  defect  structure  and  small  amount  of  polymeric  binders  and  plasticizers  [126]. 

High  explosive  can  be  detonated  by  purpose-built  detonators  [2]  or  by  high-velocity  impact  in 
accidental  scenarios  or  experimental  configuration.  The  shock  initiation  of  heterogeneous  explosive  is  a 
complicated  process  which  is  still  not  thoroughly  understood.  HMX  is  used  as  a  component  of  plastic- 
bonded  explosive  (such  as  PBX  9501),  which  is  a  typical  heterogeneous  explosive,  as  shown  in  Figure  7.1 
[127].  It  shows  PBX  9501  is  a  mixture  of  coarse  and  fine  HMX  crystal  with  diameters  between  1  and 
1000  iim ,  which  are  bound  together  by  polymers. 


Figure  7.1  Micrograph  of  PBX  9501 


When  a  shock  wave  travels  through  heterogeneous  explosive,  it  provides  heating  by  bulk  compression 
and  by  interaction  with  density  discontinuity  (such  as  voids)  and  defect  structure,  known  as  hot  spots 
[128].  On  appropriate  condition,  these  hot  spots  may  begin  to  react  and  spread  into  surrounding  explosive. 
Growing  chemical  reactions  behind  the  shock  front  raise  pressure  and  temperature  and  finally  lead  to 
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stable  detonation.  This  is  the  mesoscale  mechanism  of  detonation  of  heterogeneous  explosives  by  shock, 
as  shown  in  Figure  7.2.  According  to  Field  et  al.  [129],  the  size  of  hot  spot  ranges  from  0.1-10  /im . 


Shock  wave 


Heterogeneous  high  explosive 
Figure  7.2  Scheme  of  hot  spots  formation 

It  is  generally  accepted  that  the  initiation  of  detonation  of  heterogeneous  explosives  by  shock  can  be 
divided  into  two  distinct  stages  [128]: 

(I)  Ignition  of  a  small  fraction  of  the  explosive  at  random  sites  (zone  where  strong  shock  wave  passes 
by)  due  to  the  creation  of  hot  spots. 

(II)  Growth  to  detonation  from  the  coalescence  of  the  energy  released  from  individual  hot  spots. 

To  numerically  simulate  the  shock  initiation  of  high  explosive  in  macroscopic  scale,  appropriate 
models  need  to  be  established. 

A  classic  one-dimensional  Zeldovich-von  Neumann-Doring  (ZND)  detonation  model  is  often  used  to 
describe  the  detonation  of  heterogeneous  high  explosive  by  shock  wave  [130-132],  as  shown  in  Figure  7.3: 
first,  an  infinitely  thin  shock  wave  travelling  at  subsonic  speed  compresses  the  explosive  to  a  high 
pressure  called  von  Neumann  spike.  The  sudden  change  in  pressure  initiate  the  chemical  energy  release. 
Then,  the  energy  release  accelerates  the  flow  behind  the  spike  to  local  speed  of  sound,  which  is  known  as 
CJ  state.  After  CJ  point,  the  detonation  products  expand  backward  and  form  a  rarefaction  wave. 


CJ  point 

Detonation  wave 

- > 

von  Neumann  pike 

/ 

Gaseous  products 

Reaction  zone 

/  Inert  explosive 

One-dimensional  explosive  bar 

Figure  7.3  Scheme  of  one-dimensional  ZND  detonation  model 
The  typical  particle  velocity  histories  recorded  by  embedded  gauges  during  detonation  of 
heterogeneous  high  explosive  are  shown  in  Figure  7.4.  It  can  be  seen  from  Figure  7.4  that  shock  wave  is 
continuously  strengthened  while  propagating  in  explosive  bar  till  a  stable  detonation  wave  is  formed.  The 
distance  before  a  shock  wave  transits  to  stable  detonation  wave  is  so-called  run-to-detonation  distance. 
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Particle  velocity 


Figure  7.4  Typical  particle  velocity  histories  recorded  by  embedded  gauges 
Based  on  ZND  detonation  model,  several  numerical  models  have  been  proposed  to  describe  the 
detonation  of  high  explosive.  The  most  popular  model  is  ignition  and  growth  model,  which  is 
incorporated  in  our  SPH  code. 

Compared  with  the  standard  JWL  model  which  includes  only  a  single  EoS  for  explosive  gases, 
ignition  and  growth  model  uses  JWL  EoSs  for  both  solid  explosive  and  gaseous  products  of  explosives 
[24]: 

^  Vs  =  Ase~R +  Bse~R^s  + 


Vg  =  Age  Rlava  +  Bge  K2ava  + 


RlnVn 


C0gcVgTg 


Eqn.7.1 


V  “  “  v9 

where  p  is  pressure,  V  is  relative  volume  (V  =  V /V0  =  p0/p),  T  is  temperature,  Cv  is  average  heat 
capacity,  and  A,  B,  R±,  R2  ,  0)  are  parameters  fitted  based  on  experiment,  s  represents  solid  explosive  and 
g  represents  gaseous  products  of  explosive. 

The  speed  of  sound  is  given  by  the  approximation 


dp 

C  ~  Jdp 


Eqn.  7.2 


where  p  is  pressure  of  solid  explosive  or  gaseous  products  of  explosive,  p  is  density. 

Two  physical  assumptions  are  introduced  in  the  ignition  and  growth  model,  which  are  temperature 
equilibrium  and  pressure  equilibrium,  as  shown  as  follow  [24]: 

[Vs  =  Vg 


[Ts  =  Tg 


Eqn.  7.3 
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Temperature  equilibrium  is  often  used  as  a  standard  approach  for  closure  of  the  problem.  Temperature  is 
calculated  using  the  following  formula: 


T 

1  s 


e 

(1  -X)Cys+hCys 


Eqn.  7.3a 


where  e  is  the  specific  internal  energy  of  the  particle,  Cv  is  the  specific  heat  of  the  particle,  A  is  reaction 
progress. 

The  scheme  of  the  two  physical  assumptions  is  shown  in  Figure  7.5. 


Figure  7.5  Scheme  of  pressure  equilibrium  and  temperature  equilibrium  for  an  explosive  particle 
It  is  assumed  that  each  explosive  particle  consists  of  two  phases,  one  solid  phase  and  the  other  gas 
phase.  Since  the  mass  of  an  explosive  particle  is  fixed,  the  mass  fraction  of  solid  phase  and  gas  phase 
should  be  carefully  chosen  so  that  the  pressure  of  solid  phase  and  gas  phase  will  be  equal.  The 
temperature  of  explosive  particle  is  determined  by  the  specific  internal  energy  and  the  average  specific 
heat. 

The  reaction  rate  equation  in  ignition  and  growth  model  has  the  following  form: 

—  =  +  A2  +  A3  Eqn.  7.4 


Ui  =  I  a  -  A)*  -  1  -  a)X ,  0  <  A  <  Aigmax 

I  a2  =  G,(l  -  X)cAdVy,  0  <  A  <  AGlmax  Eqn-  7-4a 

l  A3  =  G2(l  —  A)e  Aa  pz ,  AG2min  <  A  <  1 

where  /  ,  a,  h,  x,  c,  d ,  y,  e ,  g ,  z,  G1?  G2,  A[gmax  ,  AGlmax,  AG2min  are  rate  constant  fitted  based  on 
experiment,  p  is  the  real-time  density  of  explosive  particle  during  simulation.  A±  is  ignition  item,  A2,  A3 
are  growth  items.  Different  hot  spot  models  lead  to  different  values  for  constant  x.  Constants  c,  d ,  e,  g 
are  related  to  the  choice  of  the  geometry  for  the  hot  spot  combustion  process.  Hot  spot  can  be  considered 
to  bum  outwards  from  to  void  center,  or  inward  over  the  total  grain  surface.  The  remaining  constants  /, 
G1?  G2,  z  are  found  by  fitting  simulated  pressure-time  records  to  experimentally  measured  embedded 
gauge  records. 
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In  the  implementation  of  ignition  and  growth  model,  the  state  of  explosive  can  be  described  by 
reaction  progress  A,  which  is  a  number  between  0  and  1.  0  represents  inert  explosive,  and  1  represents 
fully  reacted  explosive  products.  The  number  between  0  and  1  represents  explosion  products  in  mixing 
state. 


7.2  Molecular  dynamics  with  reactive  force  field 

Molecular  dynamics  method  is  a  powerful  tool  to  investigate  the  physical  movements  of  atoms  and 
molecules.  The  trajectories  of  atoms  and  molecules  are  determined  by  Newton’s  equations  of  motion. 

The  idea  of  molecular  dynamics  is  derived  from  quantum  chemistry  (QC),  which  is  used  to  predict 
geometries,  energies  and  vibrational  energies  for  small  molecules.  However,  QC  is  not  practical  for 
studying  the  dynamic  properties  of  larger  molecules  and  solids.  It  would  be  useful  to  have  accurate  force 
fields  to  quickly  evaluate  the  forces  and  other  dynamical  properties  such  as  the  effects  of  mechanical 
shock  waves  or  diffusion  of  small  molecules  in  polymer  and  mesoporous  zeolites  [18].  Indeed,  a  number 
of  FFs,  such  as  MM3  [133-135],  provide  accurate  predictions  of  geometries,  conformational  energy 
differences  and  heat  of  formation.  Generic  FF  such  as  DREIDING  [136]  and  UFF  [137]  allow  predictions 
for  broad  classes  of  compounds,  particularly  when  coupled  to  charge  equilibration  (QEq)  [138]  or  other 
methods  for  predicting  charges.  However,  those  force  fields  are  not  capable  to  describe  chemical 
reactivity.  Brenner  force  field  is  an  exception  because  it  is  formulated  to  describe  bond  breaking. 
However,  it  does  not  include  van  der  Waals  and  Coulomb  interactions  which  are  critical  in  predicting 
structures  and  properties  of  many  systems.  In  addition,  the  Brenner  force  field  does  not  repair 
fundamental  problems  in  the  shapes  of  dissociation  and  reactive  potential  curves.  A  method  that  has  the 
full  chemistry  of  the  breaking  and  forming  bonds  and  a  proper  description  of  the  fully  bonded  equilibrium 
geometry  of  complex  molecules  is  in  great  need. 

Thus,  a  general  bond-order-dependent  potential  named  ReaxFF  is  developed  in  which  van  der  Waals 
and  Coulomb  forces  are  included  and  the  dissociation  and  reaction  curves  are  derived  from  QC 
calculations  [18].  With  ReaxFF,  quantum  phenomena  such  as  resonance,  unsaturated  valences  in  radical 
systems  and  chemical  reactivity  can  be  accurately  described. 

Similar  to  empirical  nonreactive  force  field,  the  reactive  force  field  divides  the  system  energy  into 
various  energy  contributions  (atom  under/overcoordination,  valence  angle  terms,  torsion  angles, 
conjugated  systems,  nonbonded  van  der  Waals  interactions,  Coulomb  interactions),  as  follow  [18]: 


U  _  /7 

^cud-pm  Eh 


where  Ebond  is  the  bond  energy,  Eover  and  Eun der  is  the  energy  penalty  caused  by  over-/undercoordinated 


atom,  Fval  is  the  energy  contribution  from  valence  angle  terms,  Epen  is  the  additional  energy  penalty  for 


system  with  two  double  bonds  sharing  an  atom  in  a  valency  angle,  £j-ors  is  the  energy  of  torsion  angle, 
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fconj  is  the  energy  contributed  by  conjugation  effects,  EydWaals  is  the  energy  resulting  from  van  der 
Waals  interactions,  ^coulomb  is  the  energy  resulting  from  Coulomb  interactions. 

Since  the  introduction  of  ReaxFF,  it  has  been  successfully  applied  in  a  variety  of  reactive  dynamics 
simulations  of  hydrocarbon  organic  compounds  [139],  polymers  [21],  energetic  materials  [140-144], 
shock  process  [142,  143,  145]  and  combustion  [139,  146].  However,  when  dealing  with  solids,  ReaxFF 
does  not  account  adequately  for  long-range  London  dispersion  (van  der  Waals  attraction)  which  is 
significantly  important  in  molecular  solids,  resulting  in  equilibrium  volumes  about  10—15%  higher  than 
its  normal  density,  which  is  unacceptable.  Therefore,  an  additional  vdW-like  interaction  is  introduced  in 
the  original  ReaxFF  to  account  for  long-range  London  dispersion  [19].  The  corrected  reactive  force  field 
is  named  ReaxFF-lg. 

In  ReaxFF-lg,  the  total  energy  of  the  system  can  be  expressed  as  [19]: 

^Reax-lg  —  ^Reax  ^lg  Eqn.  7.6 

where  £Reax  is  the  energy  evaluated  from  previous  ReaxFF,  Eig  is  the  long-range-correction  terms. 

In  the  work,  an  updated  ReaxFF-lg  (ReaxFF-d3)  trained  by  QC  is  used  to  calculate  the  EoSs  of  solid 
explosive  and  gaseous  products  of  /?  —  HMX. 

7.3  Smoothed  particle  hydrodynamics  method 

7.3.1  Function  approximation  in  SPH 

SPH  was  formulated  to  solve  hydrodynamics  problems  that  are  governed  in  form  of  partial  differential 
equations  (PDEs)  with  primary  field  variables  of  density,  velocity,  and  energy  [36].  It  was  also  extended 
to  solve  problems  such  as  solid  mechanics.  In  SPH  method,  the  entire  domain  is  discretized  by  a  finite 
number  of  particles  that  carry  individual  mass  and  occupy  individual  space.  Any  field  function /(x)  in 
the  problem  domain  can  be  approximated  in  SPH  as  follow  [38]: 

(/(*))  =  f(x')W(x  -  x',  h)dx'  Eqn.  7.7 

where  H  is  problem  domain,  /  is  a  field  function  related  to  the  three-dimensional  position  vector  x, 
W(x  —  x',h )  is  smoothing  kernel  function,  h  is  smoothing  length  which  defines  influence  area  of  the 
smoothing  function  W.  Various  smoothing  kernel  functions  have  been  proposed  with  the  development  of 
SPH  method,  such  as  cubic  spline  kernel  [147],  Gaussian  kernel  [148],  quartic  and  quantic  splines  [149]. 

Smoothing  function  W  should  be  chosen  to  satisfy  a  number  of  conditions.  The  first  is  called  the 
normalization  condition: 

fn  W(x  -  x',  h)dx'  =  1  Eqn.  7.8 

The  second  is  called  the  Delta  function  property: 

lim^o  W(x  —  x',h )  =  8 (x  —  x ')  Eqn.  7.9 

The  last  is  called  the  compact  condition: 
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W (x  —  x',h )  =  0  when  \x  —  x'\  >  Kh  Eqn.  7.10 

where  k  is  a  scalar  factor  associated  with  smoothing  function  Kh  which  defines  the  effective  area  of  the 
smoothing  area.  The  effective  area  is  called  support  domain  for  the  point  at  x. 

Gaussian  kernel  is  chosen  for  the  SPH  simulation  in  the  work,  which  is  written  as: 

W(R,h)  =  ade~R2  Eqn.  7.11 

where 


R  = 


Eqn.  7.11a 


C 


Ud  =  \ 


1 

nh2  ’ 
1 


one  —  dimenional 
two  —  dimensional 
three  -  dimensional 


Eqn.  7.11b 


\n3/2h3  } 

Based  on  product  rule  for  derivatives  and  Gauss’s  Law,  the  derivative  of  function  f(x)  can  be  derived 
as  follow: 


(V  ■  /(*)>  =  Js  f(x')W  (x  -  x’,  h)  ■  ndS  -  Jn  fix')  ■  VW{x  -  x',  h)dx'  Eqn.  7.12 
where  n  is  the  unit  vector  normal  to  the  surface  S .  If  the  support  zone  is  within  problem  domain, 


according  to  compact  condition  of  smoothing  kernel  function,  we  have: 

W(x-  x',h)\s  =  0 

Eqn.  7.13 

Therefore, 

fix')  W(x  -  x' ,  h)  =  0 

Eqn.  7.14 

Thus,  Eqn.  7.12  can  be  rewritten  as: 

(V  ■  fix))  —  —  fn  fix')  ■  SIWix  —  x',  h)dx' 

Eqn.  7.15 

To  calculate  the  integration  of  Eqn.  7.7  and  Eqn.  7.15,  a  summation  over  all  particles  within  the 
support  domain  is  implemented,  as  shown  in  Figure  7.6. 


Figure  7.6  Support  zone  of  particle  i 


62 


Eqn.7.16 


Therefore,  Eqn.  7.7  and  Eqn  7.15  can  be  rewritten  as: 

_  VJV  m 


</(*;)>  = 

<v  ■  /(*;)>  =  -  Iy=i^/Oy)  ■  VWij 


Eqn.  7.17 


Eqn.  7.17a 


where  m  is  the  mass  of  particle,  p  is  the  density  of  particle,  and  Wtj  can  be  written  as 

Wtj  =  W(xi-xj,h) 

7.3.2  Discretization  of  governing  equations 

The  governing  equations  for  hydrodynamics  problems  are  Navier-Stokes  equations.  The  standard 
Lagrangian  form  of  the  Navier-Stokes  equations  consists  of,  first,  the  continuity  equation  [38]: 


Second,  the  momentum  equation: 


dp  _  dv& 

dt  ^ dx& 

Eqn.  7.18 

dva  1  doa$ 

Eqn.  7.19 

dt  p  dx& 

de  oa@  dva 

Eqn.  7.20 

dt  p  dx$ 

Third,  the  energy  equation: 


where  p  is  density,  e  is  specific  internal  energy,  va  is  velocity,  and  oa P  is  the  total  stress  tensor,  which 

are  all  dependent  variables.  xa  is  spatial  coordinate  and  t  is  time  ,  which  are  independent  variables.  The 

total  stress  tensor  oaP  consists  of  two  parts,  one  is  isotropic  pressure  p  and  the  other  is  shear  stress  xa^\ 

Gap  _  —pgaf]  _j_  Tap  Eqn.  7  21 

where  8aP  is  Kronecker  delta,  which  can  be  written  as 

'0  if  a  =£  /? 
if  a  =  /3 

For  explosive  gas,  the  viscous  shear  stress  can  be  neglected  compared  with  the  isotropic  pressure.  For 
solid  materials,  shear  stress  can  be  obtained  from  constitutive  equations  of  material.  In  implementation  of 
SPH  method,  to  approximate  density  of  particle,  an  approach  called  continuity  density  is  adopted.  After 
manipulation,  Eqn.  7.18  can  be  rewritten  as: 


8“P  = 


Eqn.  7.21a 


dp  _  w  mj  J3  dWij 

dt  Pi  Lj  =  l  p.  Vij  dxt3 

Eqn.  7.22 

can  be  respectively  rewritten  as  follow: 

d<_yJy  m.(Pi  +  Pi+n.ywU 
dt  Pj  +  Ulj)  dxf 

Eqn.  7.23 

-  1  yw  m.(Vi  a.  p>  a.  n.  VP  dWii 
-2Lj=1  mj[p2  +  p2+  Ihj)  Vij  dxp 

Eqn.  7.24 

where  Flu  is  artificial  viscosity,  /),-  is  hydrostatic  pressure  of  particle  i,  which  will  be  obtained  from  EoS. 


63 


Considering  the  release  of  chemical  heat  in  detonation,  a  source  item  is  added  to  the  energy  equation. 
Thus,  the  modified  energy  equations  used  in  SPH  can  be  expressed  as  follow: 


dfi-lvw  «,7pL.P7  .  n 


Eqn.  7.25 


where  q  is  heat  release  rate,  which  plays  critical  role  in  transition  of  shock  wave  to  detonation  wave. 
Leap-frog  algorithm  is  applied  to  integrate  the  governing  equations  over  time. 

7.3.3  Artificial  viscosity 

In  order  to  suppress  unphysical  fluctuation  in  the  simulation  of  discontinuity  such  as  shock  waves, 
artificial  viscosity  is  introduced  in  the  beginning  development  of  SPH.  Actually,  shock  wave  is  not  real 
physical  discontinuity  but  a  very  narrow  zone  where  the  status  of  medium  experiences  sudden  change.  To 
conserve  mass,  momentum  and  energy  across  the  shock  wave  front,  the  simulation  of  transformation  of 
kinetic  energy  to  heat  is  required.  Physically  the  transformation  is  done  in  the  form  of  energy  dissipation, 
which  inspires  researchers  to  introduce  artificial  viscosity  in  SPH  [150-153].  Monaghan  type  artificial 
viscosity  77^  is  implemented  in  the  following  SPH  simulation,  which  can  be  written  as: 


Vij  ■  xtj  <  0 

Vjj  '  xij  >  0 


where 


*£/  +<Pl 


cij  ~  2  ci) 

Pa  =\(pt  +  pj) 
hu  =  —  (hi  +  h: ) 


Eqn.  7.26 


Eqn.  7.26a 


Eqn.  7.26b 
Eqn.  7.26c 


Eqn.  7.26d 


Vij  =  Vt-Vj 


Eqn.  7.26e 


xij  =  xi  ~  xj  Eqn.  7.26f 

where  a  and  /?  are  constants,  <p  is  used  to  prevent  numerical  divergences  between  two  particles,  c  and  v 
represent  the  speed  of  sound  and  velocity  of  particle,  respectively.  The  viscosity  term  associated  with  a 
produces  a  bulk  viscosity,  while  the  second  term  associated  with  (>  is  intended  to  suppress  particle 
interpenetration  at  high  Mach  number.  Generally,  a,  (3,  (()  are  set  to  be  1,  1  and  0.1  respectively. 

7.3.4  Smoothing  length  update 

Smoothing  length  h  plays  significant  role  in  SPH  simulation,  which  has  direct  influence  on 
computation  efficiency  and  accuracy  of  the  solution.  If  h  is  too  small,  there  may  be  insufficient  particles 
in  the  support  domain  to  exert  force  on  a  given  particle.  Instead,  if  h  is  too  large,  all  details  of  the  particle 
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or  local  properties  may  be  smoothed  out,  and  computation  efficiency  is  reduced,  too.  The  particle 
approximation  in  SPH  depends  on  sufficient  number  of  particles  within  the  support  of  domain  of  a 
particle. 

In  the  early  development  of  SPH,  the  global  particle  smoothing  length  was  used  which  depended  on  the 
initial  average  density  of  the  system.  Later,  variable  smoothing  length  is  adopted  to  maintain  consistent 
accuracy  throughout  the  space  [148-153]. 

The  following  method  is  used  to  update  the  smoothing  length  of  each  particle  during  the 
implementation  [38]: 


d h?  _  /if  dpf 
d t  pfd  d t 


Eqn.  7.27 


where  /if,  pf are  smoothing  length,  density  of  particle  i  at  time  step  n.  For  convenience,  Eqn.  7.27  can  be 
further  written  in  the  following  form: 


d/if 


__L  —  _  _2L_  yJy  rn  .  vn  .  y  n  w  . 
d t  ~  P?dLJ=imJViJ  Vi  Wij 

where  vf  and  Vf are  velocity  and  kernel  gradient  of  particle  i  at  time  step  n. 

Thus,  the  smoothing  length  at  the  next  step  becomes 


Eqn.  7.28 


hf+1  =  hf+^At 

1  1  d  t 


Eqn.  7.29 


7.3.5  Constitutive  model  of  copper 

Johnson-Cook  (JC)  model  is  a  popular  empirical  constitutive  model  which  was  proposed  by  Johnson 
and  Cook  in  1985  [154]  to  account  for  factors  that  are  significant  when  materials  are  loaded  into  plastic 
range  in  shock  or  impact.  The  factors  include  large  strain,  large  strain  rates,  high  pressure  and  high 
temperature. 

For  materials  with  incompressible  plasticity  (such  as  copper),  the  von  Mises  stress  is  given  by: 

Gy  =  ^>]~2  Eqn.  7.30 
where  J2  stands  for  second  deviatoric  stress  invariant  that  can  be  written  as: 

h=\  [Oi  -  02)2  +  (02  -  °3)2  +  (03  -  °i)2]  Eqn-  7-30a 

where  al9  a2,  g3  are  principal  stresses. 

According  to  Johnson-Cook  model,  the  flow  stress  (yield  stress)  can  be  expressed  as  a  function  of 
plastic  strain,  strain  rate  and  temperature  as  follow: 

Gy  =  (yH-S(£p)n)(l  +  Cln£p*)(l-(r*)m)  Eqn.  7.31 

where  A,  B,  C,m  are  material  constants,  £p  is  equivalent  plastic  strain,  £p*  is  the  normalized  equivalent 
plastic  strain  rate  that  can  be  expresses  as: 
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£v*  =  —  Eqn.  7.31a 

y  £pO 

where  £p  is  equivalent  plastic  strain  rate,  £p0  is  equivalent  plastic  strain  rate  of  the  quasi-static  test.  7*  is 
expressed  as: 

Eqn.  7.31b 

Tm-To  H 

where  70  is  reference  temperature  and  Tm  is  reference  melt  temperature. 

7.3.6  Mie-Gruneisen  EoS  of  copper 

Traditionally,  an  EoS  is  equation  which  describes  the  pressure  of  a  material  as  a  function  of  density 
and  temperature.  It  is  needed  to  solve  the  continuity,  momentum  and  energy  equations  that  govern  the 
thermodynamics  transition  of  a  material  [155].  For  high-velocity  impact  problems,  thermodynamics 
states  often  change  so  rapidly  that  there  is  little  time  for  heat  transfer.  Thus  it  is  reasonable  to  consider 
high-velocity  impact  an  adiabatic  process.  Mie-Gruneisen  EoS  is  often  used  in  hydrocodes  to  simulate 
shock-compressed  solids,  which  can  be  expressed  as  [156]: 


Pofo_jP 
_  J  (l-s<p)2 


[i-±(v0-v)]  +  ge-,  ifv< 


Eqn.  7.32 


where  Q  is  Gruneisen  parameter,  c0  is  reference  sound  velocity,  p0  is  reference  density,  v  is  specific 
volume,  s  is  constant,  cp  is  expressed  as: 

(p  =  1  -  —  Eqn.  7.32a 

^0 

where 


v  =  -  Eqn.  7.32b 

P 

7.4  Benchmark  test  of  the  SPH  method  with  artificial  viscosity 

A  one-dimensional  shock  tube  problem  is  used  as  benchmark  to  validate  the  SPH  algorithm.  The 
initial  conditions  of  the  shock  tube  problem  is  described  as  [157] 

(p  =  l,v  =  0,  e  —  2.5,  p  =  1  when  x  <  0 

[p  =  0.25,  v  =  0,e  =  1.795, p  =  0.1795  whenx  >  0  Eqn.  7.33 

800  particles  are  used  in  the  simulation.  All  particles  have  the  same  mass  ofrn;  =  0.001875.  640 
particles  are  distributed  in  the  high-density  area  ranging  from  -0.6  to  0,  while  160  particles  are  distributed 
in  the  low-density  area  ranging  from  0  to  0.6.  Time  step  is  0.005s  and  total  steps  are  40.  Artificial 
viscosity  is  applied  when  using  the  traditional  SPH.  The  comparison  of  theoretical  solution  with  results 
obtained  by  traditional  SPH  is  shown  in  Figure  7.7. 
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(a)  Comparison  of  density  calculated  by  traditional  SPH  with  exact  solution 


(b)  Comparison  of  pressure  calculated  by  traditional  SPH  with  exact  solution 


67 


(c)  Comparison  of  velocity  calculated  by  traditional  SPH  with  exact  solution 


(d)  Comparison  of  specific  internal  energy  calculated  by  traditional  SPH  with  exact  solution 
Figure  7.7  One-dimensional  shock  tube  problem.  Comparison  of  characteristics  calculated  by 

traditional  SPH  with  exact  solution 


It  can  be  seen  from  Figure  7.7(a),  (b)  that  there  is  obvious  fluctuation  at  the  shock  front  when  using 
traditional  SPH,  which  indicates  that  artificial  viscosity  cannot  thoroughly  eliminate  fluctuation. 
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7.5  Detonation  of  a  ID  PBX  9501  bar 

The  following  experiment  [24]  is  studied  using  proposed  SPH  method:  A  100  mm  diameter,  12.5  mm 
thick  aluminum  flyer  plate  impacts  a  target  consisting  of:  a  90  mm  diameter,  6  mm  thick  aluminum  plate, 
a  90  mm  diameter,  20  mm  thick  PBX  9501  charge,  and  a  90  mm  diameter,  6  mm  thick  aluminum  back 
plate.  Scheme  of  the  experimental  configuration  is  shown  in  Figure  7.8.  Pressure  gauges  are  embedded  in 
the  explosive  to  record  the  variation  of  pressure  during  detonation. 

A1  fiver  plate 

^  A 1  p  late  A 1  p  late 


Figure  7.8  Scheme  of  the  configuration  of  flyer  plate  impact  experiment 

7.5.1  Case  1 

The  experiment  mentioned  above  is  simplified  to  a  one-dimensional  PBX  9501  model  as  shown  in 
Figure  7.9.  The  simplified  one-dimensional  PBX  9501  model  is  0.02  m  long.  On  the  left  there  is  a  rigid 
wall,  the  explosive  bar  will  impact  against  the  rigid  wall  at  certain  velocity.  In  all  the  cases,  particles  are 
distributed  uniformly  and  smoothing  length  of  each  particle  is  set  to  be  1.5  times  the  distance  between 
two  neighboring  particles.  The  model  consists  of  4000  particles. 


/ 

Rigid  wall  / 


PBX  9501  bar 


10  mm 


Figure  7.9  Scheme  of  simplified  flyer  plate  impact  experiment 


The  parameters  of  JWL  EoS  of  solid  PBX  9501  and  gaseous  products  of  PBX  9501,  and  JWL  EoS  of 
solid  P-HMX  and  gaseous  products  of  (B-HMX  are  listed  in  Table  6.4  and  Table  6.5.  The  parameters  for 
reaction  rate  are  as  follow: 
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Table  7.1  Reaction  rate  parameters  of  PBX  9501  [24] 


Reaction  rate  parameters  of  PBX  9501 

/ 

1.4  x  1017s-1 

a 

0.0 

b 

0.667 

X 

20.0 

c 

0.667 

d 

0.277 

y 

2.0 

e 

0.333 

9 

1.0 

z 

2.0 

Gi 

130  x  10~16Pa~y s_1 

G2 

400  x  10-16Pa-zs-1 

^-ig  max 

0.3 

^Glmax 

0.5 

^G2max 

0.5 

Figure  7.10  The  variation  of  shock  wave  on  the  explosive  bar  with  time. 
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The  vx  is  set  to  be  480  /s  .  The  resulting  impact  pressure  on  detonation  end  is  around  3.4  GPa. 
Ambient  temperature  is  25  °C.  The  variation  of  shock  wave  on  the  explosive  bar  with  time  is  shown  in 
Figure  7.10,  The  time  interval  between  two  adjacent  curves  is  0.25  fis. 

It  can  be  clearly  seen  from  Figure  7.10  that  the  shock  wave  increasingly  builds  up  while  propagating 
from  left  to  right.  When  shock  wave  arrives  at  around  0.01  m ,  the  explosive  is  detonated  and  the  shock 
wave  transforms  into  a  self-sustained  detonation  wave.  The  detonation  wave  stabilizes  at  around  0.011  m. 
This  build-up  effect  is  contributed  by  the  release  of  heat  in  detonation.  According  to  experiment  [24],  the 
experimental  run-to-detonation  distance  is  1 1  mm,  which  agrees  very  well  with  simulation. 

The  distribution  of  reaction,  vx  and  density  along  the  one-dimensional  PBX  9501  bar  at  3.5  fis  is 
shown  in  Figure  7.11. 

The  comparison  of  SPH  simulation  using  original  EoSs,  fitted  EoSs,  FVM  [24]  and  experiment  [24]  is 
shown  in  Figure  7.12,  where  the  legend  represents  the  distance  of  pressure  gauge  from  detonation  end. 

It  can  be  seen  from  Figure  7.12(b)  that  the  calculated  arrival  time  of  shock  wave  is  very  close  to 
experiment.  The  peak  pressure  calculated  from  fitted  EoS  is  around  40  GPa ,  which  also  agrees  well  with 
experiment.  However,  all  three  simulations  deviate  from  experiment  somehow.  It  probably  results  from 
three  reasons:  first,  detonation  of  high  explosive  is  very  complicated,  though  ignition  and  growth  model 
can  accurately  predict  arrival  time  and  peak  pressure  of  shock  wave,  it  is  not  capable  to  present  every 
detail  of  pressure  history;  second,  the  boundary  conditions  of  the  numerical  model  are  simplified 
compared  with  experiment;  third,  the  model  is  one-dimensional,  which  is  different  from  the  experiment  in 
reality. 


% 

cr 


(a)  Reaction 


(b)  vx 
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(c)  Density 

Figure  7.11  Reaction,  vx  and  density  of  the  one-dimensional  PBX  9501  bar  at  3.5  /is 


Time  (*ts) 

(a)  SPH  Simulation  using  original  EoSs 
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Time  (ps) 

(b)  SPH  Simulation  using  fitted  EoSs 
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(d)  Experiment 


Figure  7.12  Comparison  of  SPH  Simulation  using  original  EoSs,  fitted  EoSs,  FVM  simulation  and 

experiment. 

7.5.2  Case  2 

Configurations  of  the  simulation  are  similar  to  case  1,  except  that  the  vx  is  set  to  be  432  m/s.  Ambient 
temperature  is  50  °C.  Fitted  EoSs  obtained  from  MD  simulations  are  used.  The  evolution  of  shock  wave  is 
shown  in  Figure  7.13. 


Figure  7.13  The  variation  of  shock  wave  on  the  explosive  bar  with  time.  vx  is  set  to  be  480  /s  . 

Ambient  temperature  is  50  °C. 


74 


It  can  be  seen  from  Figure  7.13  that  due  to  the  decrease  of  impact  velocity  (thus  impact  pressure  is 
reduced),  the  shock  wave  travels  longer  before  transforming  to  detonation  wave.  According  to  experiment 
[24],  the  run-to-detonation  distance  is  13  mm,  which  agrees  well  with  simulation.  The  comparison  of 
SPH  simulation  using  original  EoSs,  fitted  EoSs,  FVM  [24]  and  experiment  [24]  is  shown  in  Figure  7.14. 


Time  his) 


(b)  SPH  simulation  using  fitted  EoSs 
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Time  fris) 


(c)  Simulation  using  FVM 


Time  (ps) 


(d)  Experiment 


Figure  7.14  Comparison  of  SPH  simulation  using  original  EoSs,  fitted  EoSs,  FVM  and  experiment 
It  can  be  seen  from  Figure  7.14  that  the  arrival  time  of  shock  wave  calculated  by  SPH  simulation  using 
fitted  EoSs  agrees  very  well  with  experiment.  However,  all  peak  pressure  (around  40  GPa )  obtained  in 
the  three  simulations,  either  SPH  with  original  EoSs,  fitted  EoS  or  FVM,  deviate  from  experiment  (30 
GPa )  a  lot.  It  indicates  that  though  ignition  and  growth  model  is  the  most  advanced  model  till  today,  it  is 
still  limited  to  depict  all  the  details  of  detonation. 
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7.5.3  Case  3 

Configurations  of  the  simulation  are  similar  to  the  other  two,  except  that  the  vx  is  set  to  be  540  m/s. 
Ambient  temperature  is  50  °C.  Fitted  EoSs  obtained  from  MD  simulations  are  used.  The  evolution  of 
shock  wave  is  shown  in  Figure  7.15,  where  the  time  interval  between  two  curves  is  0.25  /is. 

It  can  be  seen  from  Figure  7.15  that  as  impact  pressure  increases,  the  shock  wave  travels  shorter 
distance  before  transforming  to  detonation  wave. 


Figure  7.15  The  variation  of  shock  wave  on  the  explosive  bar  with  time.  vx  is  540  m/s.  Ambient 

temperature  is  50  °C. 

The  comparison  of  SPH  simulation  using  original  EoSs,  fitted  EoSs,  FVM  [24]  and  experiment  [24]  is 
shown  in  Figure  7.16. 
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Time 

(b)  SPH  simulation  using  fitted  EoSs 
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(c)  Simulation  using  FVM 


(d)  Experiment 

Figure  7.16  Comparison  of  SPH  simulation  using  original  EoSs,  fitted  EoSs,  FVM  and  experiment 
Similarly,  the  arrival  time  of  shock  wave  calculated  by  SPH  using  fitted  EoSs  agrees  well  with 
experiment,  except  peak  pressure  and  details  of  pressure  history. 
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7.6  Influence  of  various  factors  on  SPH 

The  artificial  viscosity  in  SPH  enables  it  to  suppress  oscillation  in  hydrodynamics  problems.  It  is 
necessary  to  investigate  the  influence  of  artificial  viscosity  on  simulation  results.  In  addition,  the 
influence  of  number  of  particles  on  simulation  is  also  investigated.  Original  JWL  EoSs  are  applied  in  all 
simulations  in  the  section. 

7.6.1  Influence  of  shear  viscosity  a  on  simulation 

In  this  section,  shear  viscosity  a  is  set  to  be  1.0,  2.0,  3.0  respectively.  Bulk  viscosity  /?  is  set  to  be  1.0. 
The  model  consists  of  4000  particles  in  all.  Smoothing  length  is  fixed  and  set  to  be  1.5  times  the  distance 
between  two  adjacent  particles.  The  vx  is  set  to  be  480  m/s.  Resulting  impact  pressure  on  detonation  end 
is  around  3.4  GPa.  Ambient  temperature  is  25°C. 

The  simulation  is  shown  in  Figure  7.17(a),  (b),  (c),  where  shear  viscosity  a  set  to  be  1.0,  2.0,  3.0 
respectively. 

It  can  be  seen  from  Figure  7.17  that  shear  viscosity  a  has  significant  influence  on  simulation.  As  a 
increases  ,  the  oscillation  at  the  peak  is  obviously  suppressed.  When  a  is  set  to  be  3.0,  the  curves  become 
quite  smooth  and  the  calculated  peak  pressure  is  close  to  experiment  [24].  The  calculated  run-to- 
detonation  distance  is  not  obviously  effected  by  variation  of  a. 
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(c)  a  =  3.0,  p  =  1.0 

Figure  7.17  Influence  of  shear  viscosity  a  on  simulation.  Impact  velocity  is  480  m/s.  Ambient 

temperature  is  25°C. 


7.6.2  Influence  of  bulk  viscosity  /?  on  simulation 

In  this  section,  bulk  viscosity  /?  is  set  to  be  1.0,  2.0,  3.0  respectively.  Shear  viscosity  a  is  set  to  be  1.0. 
Other  configurations  remain  the  same.  The  simulation  is  shown  in  Figure  7.18(a),  (b),  (c). 
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(a)  a  =  1.0,  p  =  1.0 


(b)  a  =  1.0,  P  =  2.0 


82 


(c)  a  =  1.0,  p  =  3.0 

Figure  7.18  Influence  of  bulk  viscosity  (3  on  SPH  simulation.  Impact  velocity  is  480  m/s.  Ambient 

temperature  is  25  °C. 

It  can  be  seen  from  Figure  7.18  that  bulk  viscosity  (3  has  much  less  influence  on  simulation  than  shear 
viscosity  a.  Though  (3  increases  from  1  to  3,  there  is  no  significant  reduction  of  oscillation  at  the  peak. 
The  calculated  run-to-detonation  distance  is  not  obviously  affected  by  variation  of  (3. 

7.6.3  Influence  of  number  of  particles  on  simulation 

In  this  section,  shear  viscosity  a  is  set  to  be  3.0,  bulk  viscosity  [3  is  set  to  be  1.0.  The  models  consist  of 
2000,  4000,  8000  particles  respectively.  The  other  configurations  are  the  same  as  mentioned  in  the 
previous  sections.  The  simulation  is  shown  in  Figure  7.19(a),  (b),  (c). 
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(a)  Number  of  particles  is  2000 


(b)  Number  of  particles  is  4000 
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(c)  Number  of  particles  is  8000 

Figure  7.19  Influence  of  number  of  particles  on  simulation.  Impact  velocity  is  480  m/s.  Ambient 

temperature  is  25  °C. 

It  can  be  seen  from  Figure  7.19  that  the  number  of  particles  has  substantial  influence  on  the 
smoothness  of  the  calculated  curves.  As  the  number  of  particles  increases,  the  calculated  curves  become 
smoother.  However,  increasing  particles  can  also  decrease  computational  efficiency.  Thus  it  needs  to  be 
balanced  between  accuracy  and  computational  efficiency.  According  to  Figure  7.19,  4000  particles  is  a 
reasonable  choice  for  the  one-dimensional  model. 

7.7  2D  QM-MD-SPH  simulation  of  high  explosives 

The  EoSs  of  solid  and  gaseous  of  products  of  HMX  obtained  from  MD  simulation  are  applied  in  a 
two-dimensional  SPH  model  of  PBX  9501  bar.  The  fitted  parameters  of  ignition  and  growth  model  are 
listed  in  Table  3.4  and  Table  3.5.  The  model  is  0.45  m  in  length,  0.025  m  in  radius.  Initial  smoothing 
length  is  set  to  be  1.5  times  the  distance  between  two  neighboring  particles.  The  vx  of  leftmost  particles 
are  set  to  be  1000  m/s  to  detonate  the  explosive  bar,  as  shown  below.  It  should  be  noted  that  EoSs  of 
solid  and  gaseous  of  products  of  HMX  obtained  from  MD  simulation  are  applied  in  the  model. 
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Figure  7.20  Two-dimensional  SPH  simulation  using  calculated  EoSs. 


The  profile  of  the  two-dimensional  SPH  model  at  different  time  is  shown  in  Figure  7.21.  The 
calculated  peak  pressure  (~43  GPa )  agrees  well  with  experiment  [24]. 
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Figure  7.21  The  profile  of  the  two-dimensional  SPH  model  at  different  time 


The  distribution  of  reaction,  vx,  density  and  smoothing  length  calculated  by  SPH  at  50  fis  is  shown  in 
Figure  7.22. 
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(c)  Density 
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Figure  7.22  The  distribution  of  reaction,  vx,  density  and  smoothing  length  at  50  fis 


7.8  QM-MD-SPH  simulation  of  aluminized  explosives 

The  addition  of  metal  particles  to  energetic  materials  is  a  well-known  means  to  improve  detonation 
performance.  Aluminum  (Al)  powders  are  widely  used  in  pyrothchnics,  rocket  propellants  and  explosives, 
as  shown  in  Figure  7.23.  Al  is  added  to  propellants  to  increase  thrust,  while  in  explosive  it  enhances  air 
blast  and  incendiary  effect. 
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Figure  7.23  Fine  aluminum  powder  (diameter  ~20  iim). 

Aluminized  explosives  are  high  explosives  using  micrometer-scale  (5~10  iim)  aluminum  particles  as 
additive.  The  mass  fraction  of  aluminum  particles  ranges  from  20%  to  40%  or  even  higher.  Aluminized 
explosives  feature  both  fast  detonation  and  slow  metal  combustion.  The  detonation  of  high  explosives 
occurs  in  microseconds,  however,  the  combustion  of  aluminum  particles  generally  takes  milliseconds  or 
even  longer  to  finish.  Thus,  combustion  of  A1  particles  in  aluminized  explosives  is  assumed  to  occur 
behind  detonation  front,  in  other  words,  during  the  expansion  of  gaseous  products.  A1  particles  act  as  inert 
ingredients  in  the  beginning  stage  of  detonation  [158].  According  to  the  experiment  conducted  by 
Trzcinski  et  al.  [159,  160],  detonation  velocity  measurement  of  aluminized  RDX  shows  that  A1  behaves  in 
the  reaction  zone  in  the  same  manner  as  a  chemical  inert  admixture.  Therefore,  different  from  ideal 
explosives  like  TNT,  HMX,  RDX,  TATB,  aluminized  explosive  is  classified  as  non-ideal  explosive.  The 
high  degree  of  secondary  exothermal  reactions  occurring  when  the  detonation  products  expand  behind  the 
detonation  zone  is  a  characteristic  feature  of  non-ideal  explosive.  Due  to  the  addition  of  aluminum 
particles,  aluminized  explosives  have  relatively  low  brisance  but  high  blast  potential  [25]. 

The  blast  performance  of  energetic  materials  can  be  measured  by  P-I  diagram,  where  P  is  peak 
pressure  and  I  is  the  impulse,  which  is  the  time  integral  of  the  pressure  generated  by  energetic  materials. 
Explosion  will  be  enhanced  by  Al  powder  since  pressure  will  be  prolonged  in  time  compared  with 
ordinary  explosive.  In  addition,  the  heated  surrounding  air  by  gaseous  products  will  also  increase  the 
pressure.  When  a  warhead  is  detonated,  for  example,  in  a  confined  room,  the  structure  will  first 
experience  a  shock  loading  then  a  quasi-static  pressure,  which  is  determining  factor  for  the  structure 
damage.  Optimal  performance  will  be  achieved  when  quasi  static  pressure  is  sufficiently  high  to  breach 
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the  structure.  When  gaseous  products  expand  and  turbulently  mix  with  surrounding  air,  the  temperature  of 
the  gases  will  drop  fast  and  the  combustion  will  be  terminated.  Thus,  refined  metal  particles  are  preferred 
because  they  burn  faster. 

The  detonation  of  aluminized  explosive  is  very  complicated  phenomenon.  The  process  includes 
reactions  between  gaseous  products  with  A1  powder,  A1  powder  with  surrounding  air,  gaseous  products 
with  surrounding  air.  The  mechanism  of  the  chemical  reactions  is  still  not  thoroughly  understood  till 
today.  Modeling  the  detonation  and  subsequent  combustion  requires  a  detailed  knowledge  of  the  chemical 
kinetics.  Abundant  research  has  been  conducted  on  the  topics,  experimentally  and  numerically.  Wilkins  et 
al.  obtained  the  Chapman- Jouguet  pressure  and  EoS  of  PBX  9404  through  experiments  [27].  Wackerle  et 
al.  investigated  the  pressure  history  of  PBX  9404  through  planar  shock  initiation  [28].  Cook  et  al. 
measured  detonation  pressure  and  detonation  velocity  of  TNT-A1  and  other  aluminized  explosives  [25]. 
Due  to  the  complexity  of  the  phenomenon,  it  takes  longer  time  for  researchers  to  understand  the 
mechanism  theoretically. 

7.8.1  Afterburning  model  for  combustion  of  aluminum  particles 

Numerical  modeling  of  metalized  explosives  poses  great  challenge  to  simulation  community  due  to 
the  nature  of  chemically  reacting  multiphase  flow.  Afterburning  model  is  a  simplified  model  to  describe 
the  combustion  of  gaseous  products,  which  is  common  in  the  simulation  of  combustion  and  detonation. 
Miller  proposed  the  idea  to  combine  ignition  and  growth  model  with  afterburning  model  [161].  According 
to  Miller,  ignition  and  growth  model  is  used  to  describe  the  detonation  of  high  explosive,  meanwhile,  the 
afterburning  model  is  used  to  describe  the  combustion  of  gaseous  products  or  aluminum  particles.  Kuhl, 
Howard  proposed  a  different  idea,  which  introduces  temperature  in  their  afterburning  model  [31]. 
Thermodynamics  code  Cheetah  is  used  to  calculate  the  constituent  of  gaseous  products  of  explosive.  In 
addition,  the  specific  internal  energy  of  components  of  explosive  gases  is  set  to  be  negative  at  the 
beginning  of  the  simulation,  which  is  incompatible  with  ignition  and  growth  model.  The  model  has  been 
used  to  study  various  explosives  such  as  SDF1  (an  aluminized  explosive,  consisting  of  45%  C4H8N808, 
35%  Al,  20%  C4H6)  ,  SDF2  (a  polyethylene-based  fuel,  consisting  of  17%  C5H8N4012,  17%  Al,  66% 
C4H6).  The  drawback  is  that  it  is  way  too  complicated  and  requires  Cheetah  to  obtain  extra  information 
(constituent  of  gaseous  products)  needed  for  the  subsequent  simulation. 

The  model  proposed  by  Miller  is  used  in  the  work.  The  EoS  for  the  gaseous  products  is  modified  to 
the  form  [32-35,  161]: 

pg  =  Age~Rl3va  +  Bge~R29v9  +  ^(c^+gQ)  Eqn.  734 

where  p  is  pressure,  V  is  volume  ratio,  and  A,  B,  R±,  R2,  co  are  parameters  fitted  from  experiment,  Q  is 
the  energy  released  in  the  combustion  of  aluminum  particles,  a  represents  the  process  of  combustion  [33]: 
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£  =  a(l -«)'«(£) 1/6  Eqn.7.35 

where  a  is  parameter  fitted  according  to  experiment,  p0  is  atmospheric  pressure,  a  is  set  to  be  1950  in  all 
simulation  of  aluminized  explosive. 

7.8.2  Numerical  model  for  simulation  of  aluminized  explosive 

The  aluminized  explosive  we  study  is  HMX-A1,  which  consists  of  69%  HMX,  15%  Al,  7.5%  GAP, 
7.5%  FomblinD  and  1%  isonate. 

The  model  is  one-dimensional,  which  is  similar  to  the  one  described  in  Figure  7.9.  There  are  totally 
2000  particles  (0.0292 m)  in  the  model.  The  initial  velocity  of  the  aluminized  explosive  bar  is  1000  m/s. 
Afterburning  model  is  combined  with  ignition  and  growth  model  for  the  simulation  of  aluminized 
explosive.  It  is  assumed  that  combustion  of  aluminum  particles  only  occurs  during  the  expansion  of 
gaseous  products. 

The  variation  of  pressure  with  time  along  the  aluminized  explosive  bar  is  shown  in  Figure  7.24.  From 
left  to  right,  there  are  10  curves  in  total  and  the  time  ranges  from  1  p.s  to  10  ps. 


Figure  7.24  Variation  of  pressure  with  time  along  the  aluminized  explosive  bar.  Each  curve  represents 

the  distribution  of  pressure  at  specific  time. 

The  growth  of  shock  wave  can  be  clearly  seen  in  Figure  7.24.  From  left  to  right,  as  shock  wave 
propagates  in  explosive,  chemical  heat  is  gradually  released  and  the  shock  wave  is  increasingly 
strengthened.  Finally,  the  shock  wave  grows  into  self-sustained  detonation  wave.  The  calculated  peak 
pressure  is  close  to  Chapman- Jouguet  pressure  obtained  in  experiment  [162]. 
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Compared  to  Figure  7.17(c),  it  can  be  seen  from  Figure  7.24  that  peak  pressure  of  aluminized 
explosive  is  reduced  from  40  GPa  to  around  25  GPa.  This  is  because  A1  powder  does  not  engage  in 
reaction  in  the  beginning  stage  of  detonation,  it  takes  longer  for  A1  powder  to  combust  than  that  needed 
for  the  detonation  of  high  explosive.  Thus,  during  detonation  of  high  explosive,  A1  powder  does  not 
contribute  to  peak  pressure,  as  shown  in  Figure  7.25.  However,  A1  powder  is  beneficial  for  the  long- 
lasting  heat  release.  The  comparison  of  reaction  progress  of  HMX  and  aluminum  at  5  p.s  is  shown  in 
Figure  7.25. 


(a)  Reaction  of  HMX  (b)  Reaction  of  aluminum  particles 

Figure  7.25  Comparison  of  reaction  progress  of  HMX  and  aluminum  particles  at  5  p.s 

It  can  be  seen  that  compared  with  high  explosive,  the  reaction  of  aluminum  particles  is  much  slower. 
The  distribution  of  pressure,  x-velocity  and  density  at  5  |is  is  shown  in  Figure  7.26. 
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Figure  7.26  Distribution  of  pressure,  x-velocity  and  density  along  the  explosive  bar 

The  pressure  history  at  different  locations  away  from  detonation  end  is  shown  in  Figure  7.27.  It  can  be 
seen  that  the  pressure  decreases  rapidly  with  distance. 


(a)  Pressure  history  at  near  field  (b)  Pressure  history  at  far  field 

Figure  7.27  Pressure  history  of  HMX-A1  at  near  field  and  far  field 

To  investigate  the  influence  of  mass  fraction  of  A1  on  performance  of  aluminized  explosive,  series  of 
simulation  is  conducted,  as  shown  in  Table  7.2. 
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Table  7.2  Aluminized  explosives  with  different  mass  fraction  of  A1 


HMX+0%A1 

5%A1 

10%  Al 

15%A1 

20%  Al 

HMX  (Kg) 

1424.9 

1363.4 

1300.2 

1234.9 

1167.7 

Al  (Kg) 

0 

86.83 

176.28 

268.47 

363.52 

Others  (Kg) 

286.4 

286.4 

286.4 

286.4 

286.4 

All  the  simulations  have  the  same  configuration  except  the  mass  fraction  of  aluminum. 


The  influence  of  mass  fraction  of  A1  on  pressure  history  at  near  field  and  far  field  is  shown  in  Figure 
7.28. 


Time  (jas)  Time  (^s) 


(a)  0.0146m  (b)  0.152  m 

Figure  7.28  Influence  of  mass  fraction  of  A1  powder  on  pressure  history  at  near  field  and  far  field 
The  variation  of  peak  pressure  with  mass  fraction  of  A1  is  shown  in  Figure  7.29. 


Al  (%)  Al  (%) 


(a)  Peak  pressure  in  near  field  (b)  Pressure  in  far  field 

Figure  7.29  Variation  of  peak  pressure  with  mass  fraction  of  aluminum 
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It  can  be  seen  from  Figure  7.29(a)  that  with  the  increase  of  Al,  the  peak  pressure  at  near  field  is 
reduced,  which  agrees  with  prior  analysis.  In  the  far  field,  the  tendency  is  reversed,  as  shown  in  Figure 
7.29(b).  As  time  passes  by,  increasing  number  of  Al  particles  start  to  combust  and  more  heat  is  released, 
enhancing  the  pressure  at  far  field.  Second  order  curve  and  first  order  curve  are  used  to  fit  the  near  field 


pressure  and  the  far  field  pressure  respectively  as  follow. 

(p  =  — 0.0133w2  +  0.0215w  +  29.8731,  0  <  w  <  20 

)  r near  ’ 

\p far  =  0.0179w  +  0.3010,  0  <  w  <  20 


Eqn.  7.36 


The  variation  of  impulse  with  mass  fraction  of  aluminum  is  shown  in  Figure  7.30.  Second  order  curve 
is  used  to  fit  the  original  curve.  According  to  the  fitted  curve,  the  optimal  performance  (maximum 
impulse)  is  achieved  when  mass  fraction  of  Al  powder  is  44.06%.  The  fitted  formula  is  as  below: 

lfar  —  —26 w2  +  2291w  +  22547, 0  <  w  <  20  Eqn.  7.37 
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Figure  7.30  Variation  of  impulse  with  mass  fraction  of  aluminum 
The  JWL  EoSs  of  solid  and  gaseous  products  of  HMX  obtained  from  molecular  dynamics  simulations 
are  applied  to  study  the  influence  of  mass  fraction  of  Al  powder  on  aluminized  explosive,  as  shown  in 
Figure  7.31.  The  configurations  of  the  model  are  similar  to  the  previous  one. 


Table  7.3  Aluminized  explosives  with  different  mass  fraction  of  aluminum 


HMX+10%A1 

11%  Al 

12%A1 

13%A1 

14%A1 

15%A1 

HMX  (Kg) 

1300.2 

1287.3 

1274.3 

1261.3 

1248.1 

1234.9 

Al  (Kg) 

176.28 

194.50 

212.82 

231.26 

249.81 

268.47 

Others  (Kg) 

286.4 

286.4 

286.4 

286.4 

286.4 

286.4 
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Figure  7.31  Influence  of  mass  fraction  of  A1  powder  on  pressure  history  at  near  field  and  far  field. 


7.9  SPH  simulation  of  cylinder  test 

The  cylinder  test  has  long  been  the  principal  method  to  obtain  EoS  parameters  for  high  explosives.  In 
this  section,  a  two-dimensional  SPH  model  of  cylinder  test  is  built.  The  internal  diameter  of  the  cylinder- 
test  device  is  2.54  cm,  the  external  diameter  is  3.06  cm  and  the  length  is  31.00  cm.  In  the  experiment, 
radial  displacement  of  copper  wall  is  recorded. 

The  two-dimensional  model  is  shown  in  Figure  7.32(a),  where  red  particles  represent  copper  and  blue 
particles  represent  PBX  9501.  The  vx  of  PBX  9501  particles  on  the  left  side  (red)  are  set  to  be  2000  m/s , 
as  shown  in  Figure  7.32(b).  In  the  simulation,  the  radial  displacement  of  copper  wall  0.2  m  away  from 

detonation  end  is  recorded  and  compared  with  experiment. 

0.2 

0.15 
0.1 
0.05 
1  0 

-0.05 
-0.1 
-0.15 

'02  o  0.1  0.2  0.3  ""  0  0.1  02  0.3 


(a) 


x  (m)  x  (m) 

Two-dimensional  SPH  model  (b)  Initial  condition 

Figure  7.32  Two-dimensional  SPH  model  of  cylinder  test 
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The  parameters  of  JWL  EoSs  of  solid  and  gaseous  products  of  PBX  9501  are  listed  in  Table  6.1  and 
Table  6.2  respectively.  The  parameters  of  reaction  rate  equation  are  listed  in  Table  7.1. 

Johnson-Cook  constitutive  model  is  used  to  describe  the  flow  behavior  of  copper.  Johnson-Cook 
model  is  an  empirical  model  used  to  calculate  flow  stress  of  solid  materials  subject  to  large  strains,  high 
strain  rates  and  high  temperature  [154].  The  parameters  in  JC  model  for  copper  are  listed  as  below: 


Table  7.4  Parameters  in  the  JC  model  for  copper  [154] 


A(MPa) 

B(MPd) 

C 

n 

m 

TV  oom(V) 

Tme  ItOO 

89.63 

291.6 

0.025 

0.31 

1.09 

294 

1356 

In  the  implementation  of  SPH  code,  Mie-Gruneisen  EoS  is  used  to  calculate  the  pressure  of  copper 
under  shock  loading.  The  parameters  in  Mie-Gruneisen  EoS  for  copper  are  listed  in  Table  7.5  [163,  164]: 


Table  7.5  Parameters  in  the  MG  EoS  for  copper  [163,  164] 


Po(.kg/m3) 

c(m/s ) 

5 

G(GPa) 

Y 

8940 

3940 

1.49 

47.7 

1.96 

The  simulation  result  is  shown  in  Figure  7.33,  where  red  particles  represent  copper  and  blue  particles 
represent  PBX  9501.  Shear  viscosity  a  is  set  to  be  1,  bulk  viscosity  /?  is  set  to  be  1,  (p  is  set  to  be  0.1. 
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Figure  7.33  Expansion  of  gaseous  products  and  copper  wall  at  different  time 
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From  Figure  7.33  it  can  be  observed  that  with  the  expansion  of  explosive  products,  the  copper  wall  is 
gradually  dilated  and  breached. 

The  distribution  of  vx ,  density  p,  pressure  p  at  40  ps  are  shown  in  Figure  7.34. 
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Figure  7.34  Distribution  of  various  properties  at  40  ps 


The  radial  displacement  of  copper  wall  located  at  0.2  m  away  from  detonation  point  is  compared  with 
experiment  [165],  as  shown  in  Figure  7.35. 


Figure  7.35  Comparison  of  calculated  radial  displacement  of  copper  wall  with  experiment 


According  to  Figure  7.35,  the  radial  displacement  of  copper  wall  obtained  from  simulation  does  not 
agree  very  well  with  experiment.  This  is  probably  because  the  SPH  model  is  two-dimensional.  In  addition, 
accurate  simulation  of  cylinder  test  requires  more  complex  models,  such  as  failure  model  for  copper  wall. 
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Chapter  8-SPH  using  Godunov  scheme 

In  the  beginning  development  of  SPH  method,  artificial  viscosity  is  introduced  to  suppress  unphysical 
fluctuation.  Afterwards,  it  is  widely  accepted  as  a  standard  means  to  simulate  discontinuous  phenomenon 
such  as  shock  waves.  Abundant  literature  can  be  found  on  the  application  of  SPH  with  artificial  viscosity 
[166-170].  However,  the  parameters  (shear  viscosity  and  bulk  viscosity)  of  artificial  viscosity  need  to  be 
tuned  before  each  simulation,  which  can  be  quite  time-consuming.  Furthermore,  artificial  viscosity  only 
suppresses  unphysical  fluctuation  instead  of  eliminating  them  thoroughly.  Under  certain  circumstances, 
artificial  viscosity  may  not  function  and  cause  catastrophic  failure.  Therefore,  an  innovative  technique  is 
required  to  eliminate  artificial  viscosity  in  traditional  SPH  algorithm. 

The  Riemann  solver  in  computational  fluid  dynamics  (CFD)  has  inspired  researchers  a  lot.  Godunov 
proposed  a  numerical  scheme  to  solve  exact  or  approximate  Riemann  problem  in  CFD  in  1959 
[171]. Various  Riemann  solvers  have  been  proposed  since  then,  such  as  Roe  solver,  HLLE  solver,  HLLC 
solver  and  so  on  [172],  which  inspires  researchers  to  introduce  Godunov’s  method  in  traditional  SPH 
method  to  replace  artificial  viscosity.  For  example,  Inutsuka  proposed  the  idea  to  apply  Godunov’s 
method  in  SPH  method  [173,  174].  Afterwards,  Monaghan  discussed  the  feasibility  of  applying  Riemann 
solver  in  SPH  method  [175].  Inspired  by  Monaghan’s  idea,  Parshikov  rearranged  the  governing  equations 
and  introduced  approximate  Riemann  solver  in  traditional  SPH  method  and  his  method  yields  reasonably 
good  results  [176,  177].  Cha  et  al.  implemented  and  reviewed  several  types  of  Godunov  SPH  [178].  Since 
then,  Godunov  SPH  has  been  increasingly  applied  in  hydrodynamics  and  other  problems,  which  can  be 
found  in  [179-181]. 

Detonation  of  high  explosives  (such  as  HMX,  PBX  9404,  PBX  9501)  is  an  extreme  phenomenon 
containing  complicated  chemical  reactions  and  strong  discontinuity  such  as  shock  wave  [2].  It  is  not  only 
hazardous  but  also  expensive  to  conduct  experiment  on  explosives.  Thus,  numerical  investigation  of  high 
explosive  has  long  been  the  interest  of  researchers.  Tarver  and  McGuire  proposed  ignition  and  growth 
model  and  applied  the  model  in  the  study  of  detonation  waves  of  high  explosive  such  as  TNT,  TATB, 
LX-17,  PBX  9404  [3,  4,  182-185].  Since  the  proposal  of  ignition  and  growth  model,  it  has  been  the  most 
popular  mathematical  model  to  describe  detonation  of  high  explosives  and  incorporated  in  commercial 
software  such  as  ANSYS  LS-DYNA.  Kapila  et  al.  described  the  ignition  and  growth  model  in  a  finite- 
volume  method  and  applied  the  model  in  simulation  of  LX-17  [186].  Souers  described  a  simplified 
reactive  flow  model  (JWL++  model)  and  used  the  model  to  study  size  effect  of  ANFO  K1  [6].  However, 
JWL++  model  uses  simple  mixture  rule  instead  of  making  real  physical  assumptions  in  the  description  of 
detonation,  which  makes  it  much  less  accurate  than  ignition  and  growth  model.  Garcia  and  Tarver  further 
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described  a  three-dimensional  ignition  and  growth  reactive  flow  model  of  PBX  9502  and  compared  the 
simulation  with  experiment  [187].  In  this  section,  a  Godunov  SPH  which  integrates  ignition  and  growth 
model  is  proposed  to  investigate  the  detonation  of  PBX  9501.  The  proposed  method  can  be  extended  to 
other  high  explosives  such  as  PBX  9404  and  so  on. 

8.1  Governing  equations  of  Godunov  SPH 

The  governing  equations  in  traditional  SPH  (viscous  shear  stress  is  neglected  in  the  equations)  can  be 
written  as 


d p.  A  rn.p. 

ur  j=i  p. 

dv,_  fn>,(p,  +  p,) 

it  ji  P,P,  ‘  “ 

d  <?,.  1  ^mJ(pi  +  pj)  . 


Eqn.  8.1 

Eqn.  8.2 
Eqn.  8.3 


dt  2  M  pjPi 

where  Vk.  is  kernel  function,  p  is  density,  v  is  velocity,  p  is  pressure,  e  is  specific  internal  energy,  q  is 
heat  release  rate,  t  is  time. 

In  SPH,  the  kernel  H'; is  a  function  of  |/'  - r  |  jh (h  is  smoothing  length),  its  gradient  is  written  as 

•  r  ri~ri 


f\  \\ 

rr, 

h 

y 


=  W' 


Eqn.  8.4 


where 


v  = —  Eqn.  8.4a 

dt 

To  introduce  Riemann  solver  in  SPH,  it  is  assumed  that  a  discontinuity  exists  in  the  middle  of  each 
particle  pair,  for  example,  particle  i  and  particle  j  ,  as  shown  in  Figure  8.1. 
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The  velocity  and  pressure  at  the  discontinuity  can  be  computed  using  the  following  approximations 
[177]: 


R  _VjPjCj+V!<P.C,-Pj+P, 


vij  = 


pfj +Pici 


PiPfi  +  PiPPi  -pPiPPiip  -V*) 


J  J  ' J  J'  11  V  J 

PjCj+PPi 


where  c  is  sound  velocity,  v  is  written  as 


r  —r 
vR  =  v  •  — — - 

rj  ~r 


The  Riemann  solver  is  introduced  in  traditional  SPH  via  the  substitution: 


Eqn.  8.5 


Eqn.  8.6 


Eqn.  8.6a 


A  /  R  .  R\  . 

2(v<  +vJ^vu 


^P'+pj^p; 


Eqn.  8.7 
Eqn.  8.8 


Thus,  after  rearrangements,  the  continuity,  momentum,  modified  energy  equations  of  Godunov  SPH 
can  be  described  as  below  [188]: 


Eqn.  8.9 

dv,  mjPlw,rj~ri 

d t  mPjPP  y 

Eqn.  8.10 

^'=-2  ±mjP;M-v;Rw;+q 
d t  77  PjPih 

Eqn.  8.11 

8.2  Artificial  viscosity  in  traditional  SPH 

In  traditional  SPH,  to  simulate  the  transformation  of  kinetic  energy  to  heat  energy,  Monaghan  type 
artificial  viscosity  77.  is  introduced  [150,  151],  which  can  be  expressed  as 

l-aTjh  +  M  _ 


-ijrij  mj 

Pu 


V*«<0 


Eqn.  8.12 


0  ,  v ..  •  xtj  >  0 
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where 


A  -  V*** 

Eqn.  8.12a 

rij  2  2 

xa  +(P 

^  1 

II 

+ 

Eqn.  8.12b 

pij^pi+pj) 

Eqn.  8.12c 

Eqn.  8.12d 

vn=vi-vj 

Eqn.  8.12e 

XiJ=Xi-Xj 

Eqn.  8.12f 

where  a  and  /?  are  constants.  #?is  used  to  prevent  numerical  divergences  between  two  particles,  c  and 
v  represent  the  speed  of  sound  and  particle  velocity,  respectively.  The  a  and  /?  need  to  be  tuned  before 
each  simulation.  In  general  cases,  a  and  /?  are  set  to  be  1,  cp  is  set  to  be  0.1. 

8.3  Ignition  and  growth  model 

The  ignition  and  growth  model  is  integrated  in  our  SPH  code.  Compared  with  traditional  JWL  model 
which  includes  only  a  single  equation  of  state  for  explosive  gases,  ignition  and  growth  model  employs 
JWL  EoSs  for  both  solid  explosive  and  gaseous  products  of  explosives  [24]: 


pc  =  Ae 

I  s  s 


p ,  =  V" 


-R,.V,  ,  r>  V',  <O.Cv.T. 


+  Be 


Ruvg  +B  e~R2 gVg  , 


K 

(o„CV!T„ 

~V~ 


Eqn.  8.13 


where  p  is  pressure,  V  is  relative  volume,  T  is  temperature,  C\  is  average  heat  capacity,  A,  B,  R{  ,  R2 
and  co  are  parameters  fitted  based  on  experiment,  s  represents  solid  explosive  and  g  represents  gaseous 
products  of  explosive.  The  speed  of  sound  is  given  by 

\dP  Eqn.  8.14 


c  = 


I  dp 


where  p  is  pressure,  p  is  density. 

Two  physical  assumptions  are  introduced  in  ignition  and  growth  model,  which  are  temperature  and 
pressure  equilibrium,  as  shown  below 


Ps=Pg 
T  =T 


Eqn.  8.15 


Temperature  equilibrium  is  often  used  as  a  standard  approach  for  closure  in  hydrocodes. 
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The  reaction  rate  equation  is  written  as  [24] 


—  =  /(  \-X)b  — -1-a  +  G1(l-A)cAdpy+G2(l-A)eAgpz 

d  t  pn  J  - - - '  ' - - - - 


Eqn.  8.16 


are  rate  constants  fitted  based  on 

LrZmin 


experiment,  p  is  real-time  density  of  explosive  particle  during  simulation. 

In  the  implementation  of  ignition  and  growth  model,  the  state  of  explosive  can  be  described  by 
reaction  progress  X  ,  which  is  a  number  between  0  and  1 .  0  represents  solid  explosive,  while  1  represents 


fully  reacted  explosive  products.  The  number  between  0  and  1  represents  explosion  products  in  mixing 
state. 


8.4  Benchmark  test  for  Godunov  SPH 

A  one-dimensional  shock  tube  problem  is  used  as  benchmark  to  validate  the  Godunov  SPH  algorithm. 
The  initial  conditions  of  the  shock  tube  problem  is  described  as  [157] 


Eqn.  8.17 


800  particles  are  used  in  the  simulation.  All  particles  have  the  same  mass  of  —  0.001875.  640 
particles  are  distributed  in  the  high-density  area  ranging  from  -0.6  to  0,  while  160  particles  are  distributed 
in  the  low-density  area  ranging  from  0  to  0.6.  Time  step  is  0.005  s  and  total  steps  are  40.  Artificial 
viscosity  is  applied  when  using  the  traditional  SPH.  The  comparison  of  results  obtained  by  traditional 
SPH  with  that  obtained  by  Godunov  SPH  is  shown  below 


1.2 


— Exact  solution 
-©-Traditional  SPH 
-*-Godunov  SPH 


°-(j.4  -0.3  -0.2  -0.1  0  0.1  0.2  0.3  0.4 

x(m> 

(a)  Comparison  of  density  calculated  by  traditional  SPH  with  artificial  viscosity,  Godunov  SPH  with 

exact  solution 
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(b)  Comparison  of  pressure  calculated  by  traditional  SPH  with  artificial  viscosity,  Godunov  SPH  with 

exact  solution 


(c)  Comparison  of  velocity  calculated  by  traditional  SPH  with  artificial  viscosity,  Godunov  SPH  with 

exact  solution 
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3 


— Exact  solution 
-©-Traditional  SPH 
-*-Godunov  SPH 


x(m) 

(d)  Comparison  of  specific  internal  energy  calculated  by  traditional  SPH  with  artificial  viscosity, 

Godunov  SPH  with  exact  solution 

Figure  8.2  One-dimensional  shock  tube  problem.  Comparison  of  properties  calculated  by  traditional 
SPH  with  artificial  viscosity,  Godunov  SPH  with  exact  solution 

It  can  be  seen  from  Figure  8.2(a),  (b)  that  there  is  violent  fluctuation  at  the  discontinuity  when  using 
traditional  SPH.  Even  artificial  viscosity  cannot  thoroughly  eliminate  the  phenomenon.  However,  when 
using  Godunov  scheme,  the  fluctuation  is  obviously  smoothed.  It  indicates  that  Godunov  scheme  can 
effectively  eliminate  fluctuation  at  discontinuity  while  maintaining  the  accuracy  of  traditional  SPH. 
Another  advantage  of  Godunov  SPH  is  that  the  time-consuming  process  of  searching  for  appropriate 
values  for  artificial  viscosity  is  no  longer  needed. 

8.5  Detonation  of  a  ID  PBX  9501  bar 

The  following  experiment  [24]  is  studied  using  SPH  simulations: 

A  100  mm  diameter,  12.5  mm  thick  aluminum  flyer  plate  impacts  a  target  consisting  of:  a  90  mm 
diameter,  6  mm  thick  aluminum  plate,  a  90  mm  diameter,  20  mm  thick  PBX  9501  charge,  and  a  90  mm 
diameter,  6  mm  thick  aluminum  back  plate.  Pressure  gauges  are  embedded  in  the  explosive  to  record  the 
variation  of  pressure  during  detonation. 

8.5.1  Case  1 

The  simplified  one-dimensional  PBX  9501  model  is  0.02  m  long,  as  shown  in  Figure  7.9.  On  the  left 
there  is  a  rigid  wall,  the  explosive  bar  impacts  the  rigid  wall  at  certain  velocity.  Particles  are  uniformly 
distributed  and  smoothing  length  of  each  particle  is  set  to  be  1.5  times  the  distance  between  two  adjacent 
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particles.  The  model  consists  of  4000  particles.  The  distribution  of  pressure  along  the  explosive  bar  at 
0.01  iis  is  shown  in  Figure  8.3. 


20  | — - - 1 - - [— -  |  - i  |  i 


-5  - 

_A  Q  ‘  -  ‘  ‘  I  -  ‘  -  ‘  *  ‘  -  ‘  ‘  I  ‘  ‘  -  ‘  1  ‘  ■  ‘  ‘ 

0  0.005  0.01  0.015  0.02  0.025 

x(m) 

Figure  8.3  Distribution  of  pressure  at  0.01  /rs 

The  parameters  of  JWL  EoS  of  solid  PBX  9501  are  listed  in  Table  6.1.  The  parameters  of  JWL  EoS  of 
gases  products  of  PBX  9501  are  listed  in  Table  6.2.  The  parameters  of  reaction  rate  are  listed  in  Table  7.1. 

The  vx  is  set  to  be  480  m/s.  The  variation  of  shock  wave  along  the  explosive  bar  with  time  obtained 
from  Godunov  SPH  is  shown  in  Figure  8.4,  where  each  curve  represents  the  profile  of  pressure  on  the 
explosive  bar  at  a  specific  time: 


Figure  8.4  Calculated  variation  of  shock  wave  along  the  explosive  bar  with  time  by  Godunov  SPH.  vx 
is  set  to  be  480  m/s.  Ambient  temperature  is  25  °C. 
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It  can  be  clearly  seen  from  Figure  8.4  that  the  shock  wave  builds  up  while  propagating  from  left  to  the 
right.  When  shock  wave  arrives  at  around  0.011  m ,  explosive  bar  is  detonated  and  the  shock  wave 
transforms  to  self-sustained  detonation  wave.  The  detonation  wave  stabilizes  at  around  0.013  m. 
According  to  experiment  [24],  the  experimental  run-to-detonation  distance  is  11  mm  ,  which  agrees  very 
well  with  simulation. 

The  profile  of  reaction  and  density  of  the  one-dimensional  PBX  9501  bar  at  3.5  [is  is  shown  in  Figure 
8.5. 


x  (m)  x  (m) 

(a)  Reaction  (b)  Density 

Figure  8.5  Profile  of  reaction  and  density  along  the  one-dimensional  PBX  9501  bar  at  3.5  [is 

The  comparison  of  pressure  history  calculated  by  traditional  SPH,  Godunov  SPH,  FVM  [24]  with 
experiment  [24]  is  shown  in  Figure  8.6. 
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Time  (jxs)  Time  (p.s) 

(c)  FVM  (d)  Experiment 

Figure  8.6  Case  1.  Comparison  of  pressure  history  calculated  by  traditional  SPH,  Godunov  SPH,  FVM 
with  experiment.  Ambient  temperature  is  25  °C.  vx  is  set  to  be  480  m/s. 


According  to  Figure  8.6(a),  (b),  the  calculated  pressure  histories  using  traditional  SPH  and  Godunov 
SPH  are  very  similar.  Meanwhile,  the  curves  calculated  by  both  SPH  methods  are  much  smoother  than 
those  calculated  by  FVM.  It  is  because  SPH  is  Fagrangian  method,  which  means  that  it  is  easier  to  track 
the  interface  of  materials.  The  peak  pressure  calculated  by  both  traditional  SPH  and  Godunov  SPH  are 
very  close  to  40  GPa ,  which  agrees  well  with  experiment  in  Figure  8.6(d).  In  addition,  it  can  be  seen  from 
Figure  8.6(a),  (b)  that  the  run-to-detonation  time  is  around  3.5  /rs  ,  which  is  very  close  to  experiment  (3.4 
It  also  should  be  noted  that  the  results  calculated  by  traditional  SPH  and  Godunov  SPH  are  kind  of 
different  from  experiment,  which  are  probably  because:  firstly,  the  boundary  condition  of  SPH 
simulations  is  simplified  to  rigid  wall,  which  is  different  from  experimental  configuration;  secondly,  the 
SPH  model  is  one-dimension,  which  is  different  from  experiment;  thirdly,  though  ignition  and  growth 
model  is  significantly  improved  compared  with  JWF++  model  and  JWF  model,  it  cannot  present  all  the 
details  of  detonation. 

The  calculated  history  of  particle  velocity  is  shown  as  below 
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Figure  8.7  History  of  particle  velocity  calculated  by  Godunov  SPH 

It  should  be  noted  that  in  Figure  8.7  the  benchmark  of  particle  velocity  is  set  to  be  0. 

8.5.2  Case  2 

In  case  2,  all  configurations  remain  unchanged  except  that  vx  is  set  to  be  432  m/s.  The  variation  of 
shock  wave  along  the  explosive  bar  with  time  obtained  from  Godunov  SPH  is  shown  in  Figure  8.8. 


Figure  8.8  The  variation  of  shock  wave  along  the  explosive  bar  with  time  obtained  from  Godunov 
SPH.  vx  is  set  to  be  432  m/s.  Ambient  temperature  is  50  °C. 
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It  can  be  seen  in  Figure  8.8  that  due  to  the  decrease  of  impact  velocity,  the  shock  wave  travels  longer 
(at  around  0.013  m)  before  transforming  to  detonation  wave.  According  to  experiment,  the  run-to- 
detonation  distance  is  0.013  m.  The  simulation  agrees  very  well  with  experiment. 

The  comparison  of  pressure  history  calculated  by  traditional  SPH,  Godunov  SPH,  FVM  [24]  with 
experiment  [24]  is  shown  in  Figure  8.9. 


Time  (^s)  Time  (*is) 


(c)  FVM  (d)  Experiment 

Figure  8.9  Case  2.  Comparison  of  pressure  history  calculated  by  traditional  SPH,  Godunov  SPH,  FVM 
with  experiment.  Ambient  temperature  is  50  °C.  vx  is  set  to  be  432  m/s. 

It  can  be  seen  from  Figure  8.9(b),  (d)  that  the  pressure  histories  calculated  by  Godunov  SPH  are  quite 
close  to  experiment  in  some  ways,  such  as  arrival  time  of  shock  wave  and  certain  details.  According  to 
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Figure  8.9(b),  the  calculated  run-to-detonation  time  is  around  4.3  /rs,  which  agrees  well  with  experiment 
(4.3  iis). 

The  calculated  history  of  particle  velocity  is  shown  as  below 


Figure  8.10  Calculated  history  of  particle  velocity  calculated  by  Godunov  SPH 

8.5.3  Case  3 

The  configurations  are  similar  to  case  1  and  case  2  except  that  the  vx  is  set  to  be  540  m/s.  The 
calculated  variation  of  shock  wave  along  the  explosive  bar  with  time  by  Godunov  SPH  is  shown  in  Figure 


Figure  8.1 1  The  variation  of  shock  wave  along  the  explosive  bar  with  time  obtained  from  Godunov 
SPH.  vx  is  set  to  be  540  m/s.  Ambient  temperature  is  50  °C. 
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It  can  be  seen  from  Figure  8.11  that  as  impact  velocity  increases,  the  shock  wave  travels  shorter  before 
transforming  to  detonation  wave.  The  calculated  run-to-detonation  distance  is  around  0.009  m ,  which  is 
very  close  to  experiment  (0.008  m). 

The  comparison  of  pressure  history  calculated  by  traditional  SPH,  Godunov  SPH,  FVM  [24]  with 
experiment  [24]  is  shown  in  Figure  8.12. 


Time  (^s)  Time  (f.is) 


(a)  Traditional  SPH 


(b)  Godunov  SPH 


Figure  8.12  Case  3.  Comparison  of  pressure  history  calculated  by  traditional  SPH,  Godunov  SPH, 
FVM  with  experiment.  Ambient  temperature  is  50  °C.vx  is  set  to  be  540  m/s. 
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It  can  be  seen  from  Figure  8.12(b)  that  the  calculated  peak  pressure  of  gauge  located  at  10  mm  is 
around  24  GPa ,  which  is  quite  close  to  peak  pressure  recorded  in  experiment  (around  27  GPa ).  According 
to  Figure  8.12(b),  the  calculated  run-to-detonation  time  is  around  2.5  / is ,  which  is  close  to  experiment 
(2.4  fis). 

The  calculated  history  of  particle  velocity  is  shown  as  below 


8.6  2D  simulation  of  PBX  9501  using  Godunov  SPH 

A  two-dimensional  model  of  PBX  9501  bar  is  developed  to  investigate  the  stability  and  accuracy  of 
the  Godunov  scheme.  The  parameters  of  ignition  and  growth  model  are  the  same  as  those  mentioned  in 
previous  section.  The  model  is  0.3  m  in  length,  0.0127  m  in  radius.  Smoothing  length  is  set  to  be  1.5 
times  the  distance  between  two  neighboring  particles.  The  vx  of  leftmost  particles  are  set  to  be  2000  m/s 
to  detonate  the  explosive  bar,  as  shown  below 
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x{m) 


Figure  8.14  Profile  of  two-dimensional  SPH  model 


The  comparison  of  profiles  calculated  by  traditional  SPH  (left)  and  Godunov  SPH  (right)  at  different 
time  is  shown  in  Figure  8.15. 
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(a)  Comparison  of  profiles  calculated  by  traditional  SPH  (left)  and  Godunov  SPH  (right)  at  10  /is 
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(b)  Comparison  of  profiles  calculated  by  traditional  SPH  (left)  and  Godunov  SPH  (right)  at  20  fis 
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(c)  Comparison  of  profiles  calculated  by  traditional  SPH  (left)  and  Godunov  SPH  (right)  at  30  /is 
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(d)  Comparison  of  profiles  calculated  by  traditional  SPH  (left)  and  Godunov  SPH  (right)  at  40  / is 
Figure  8.15  Comparison  of  profiles  calculated  by  traditional  SPH  and  Godunov  SPH  at  different  time 


The  comparison  of  reaction,  vx  and  density  calculated  by  traditional  SPH  and  Godunov  SPH  at  40  /rs 
is  shown  in  Figure  8.16. 
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(c)  Comparison  of  vx  calculated  by  traditional  SPH  (left)  and  Godunov  SPH  (right)  at  40  ps 
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(c)  Comparison  of  p  calculated  by  traditional  SPH  (left)  and  Godunov  SPH  (right)  at  40  ps 
Figure  8.16  Comparison  of  reaction,  vx  and  density  calculated  by  traditional  SPH  and  Godunov  SPH 

at  40  ps 
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Conclusion 

Equation  of  state  (EoS)  plays  a  fundamental  role  in  numerical  modeling  of  high  explosives  based  on 
hydrodynamic  formulations.  Traditionally,  the  EoS  of  explosives  can  only  be  obtained  through 
experiments,  which  is  not  only  costly  and  time-consuming,  but  also  hazardous.  This  report  presents  a 
method  to  calculate  the  EoS  of  solid  p-HMX  and  its  gaseous  products  using  a  series  of  molecular 
dynamics  (MD)  simulations.  The  reactive  force  fields  ReaxFF-d3  used  in  these  MD  simulations  are 
obtained  through  quantum  mechanics  (QM)  simulations.  A  QM-MD-SPH  method  has  been  successfully 
developed,  coded  and  tested  using  a  number  of  benchmark  problems  and  experimental  results  for  various 
types  of  high  explosives. 

The  geometric  structure  of  p-HMX  is  obtained  from  DFT-D2  study.  It  is  found  that  the  calculated 
lattice  constants  of  p-HMX  agree  well  with  experiments.  In  addition,  the  mechanical  properties,  EoSs  of 
Aluminum  and  its  oxide  are  also  studied  through  a  series  of  first-principle  MQ  simulations.  The  stability, 
mechanical  properties  and  EoSs  of  iron  are  obtained  from  the  first  principle. 

The  calculated  EoSs  are  incorporated  in  our  in-house  QM-MD-SPH  codes.  The  simulation  using  fitted 
EoSs  is  compared  with  experimental  results  (from  the  open  literature),  and  results  obtained  by  the  other 
numerical  methods.  It  is  found  that  the  QM-MD-SPH  using  the  fitted  EoSs  from  MD  simulation  results 
can  well  predict  arrival  time  of  shock  waves  and  the  peak  pressure  in  many  cases.  However,  due  to  the 
complexity  of  detonation  problems,  all  numerical  models  can  not  accurately  predict  detailed  pressure 
histories,  though  they  can  roughly  predict  the  build-up  process.  Other  reasons  for  the  error  may  include  1) 
simplified  boundary  conditions  and  settings  in  the  numerical  model,  which  may  be  quite  different  from 
the  experimental  setup;  2)  the  purity  of  the  materials  may  also  contributes  to  the  errors;  3)  experimental 
error. 

In  addition,  the  influence  of  shear  viscosity  a  ,  bulk  viscosity ) 3 ,  and  number  of  particles  on  simulation 
in  QM-MD-SPH  is  analyzed,  through  the  simulation  of  the  explosion  of  one-dimensional  PBX  9501  bar 
model  ignited  by  impact.  It  is  found  that  with  the  increase  of  a  ,  the  fluctuation  at  peak  pressure  is 
significantly  reduced,  and  3  is  a  reasonable  choice  for  simulations.  It  is  also  found  that  bulk  viscosity  /? 
has  much  less  influence  in  reducing  fluctuation  at  peak  pressure.  Both  a  and  /?  do  not  have  obvious 
influence  on  arrival  time  of  shock  waves.  Moreover,  when  number  of  particles  increases,  the  calculated 
curves  becomes  smoother  and  more  details  of  the  explosion  events  can  be  revealed.  This  is,  of  course,  at 
the  expenses  of  the  computational  resources.  For  the  proposed  one-dimensional  PBX  9501  model,  4000 
particles  (0.02  m)  is  a  reasonable  choice  high  explosives  we  have  studied  (on  a  normal  PC  without 
parallelization).  For  the  high  resolution  solutions  of  2D  and  3D  problems,  more  powerful  computer  with 
parallelization  capabilities  is  needed. 
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In  the  study  of  aluminized  explosives,  an  afterburning  model  is  integrated  with  ignition  and  growth 
model  and  implemented  in  our  SPH  code  for  simulating  the  explosion  of  the  aluminized  explosive  HMX- 
A1  with  various  A1  mass  fractions.  It  is  found  that  in  the  near  field,  with  the  increase  of  mass  fraction  of 
aluminum  powder,  the  peak  pressure  is  reduced.  This  is  because  the  burning  of  A1  is  much  slower 
compared  with  the  detonation  of  the  explosives.  Therefore,  at  the  early  stage,  the  aluminum  powder  plays 
only  a  role  of  diluting  the  explosion  power.  In  the  far  field,  however,  the  pressure  is  strengthened.  This  is 
because  the  burning  of  the  aluminum  powder  adds  on  the  explosion  power  at  the  later  time. 

A  Godunov  SPH  method  is  incorporated  into  our  method,  using  the  ignition  and  growth  model. 
Godunov  approach  eliminates  the  need  of  using  the  artificial  viscosity  in  traditional  SPH  to  reduce  the 
oscillations  near  the  shock  front.  To  validate  the  Godunov  scheme,  a  one-dimensional  shock  tube  problem 
is  investigated  using  both  traditional  SPH  with  artificial  viscosity  and  Godunov  SPH.  It  is  found  that 
Godunov  SPH  can  effectively  reduce  unphysical  fluctuation,  without  the  use  of  artificial  viscosity.  It  is 
recommended  for  the  simulation  of  high  explosives  where  shocks  always  exist. 

The  Godunov  SPH  is  also  applied  to  investigate  the  detonation  of  a  one-dimensional  PBX  9501  bar. 
We  found  that  the  Godunov  SPH  performs  better  than  traditional  SPH  with  better  accuracy.  In  addition, 
due  to  its  advantage  in  describing  interfaces  of  materials,  the  calculated  pressure  histories  by  Godunov 
SPH  are  smoother.  Further,  a  two-dimensional  model  of  PBX  9501  bar  is  developed  to  examine  the 
accuracy  of  Godunov  SPH  with  respect  to  traditional  SPH.  It  is  found  that  some  unphysical  behavior 
(observed  on  the  outer  particles  of  the  two-dimensional  SPH  model)  in  the  traditional  SPH  simulation  is 
removed  when  using  Godunov  approach.  Overall,  Godunov  SPH  is  a  promising  method  compared  with 
traditional  SPH.  The  proposed  Godunov  SPH  method  can  be  extended  to  study  other  similar  explosives. 
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Appendix 

The  ReaxFF-d3  reactive  force  field  obtained  from  quantum  mechanics  simulation  is  as  follow: 
Reactive  MD-force  field 

39  !  Number  of  general  parameters 


50.0000 

! Comment 

here 

9.4514 

! Comment 

here 

30.3695 

! Comment 

here 

216.7028 

! Comment 

here 

12 .7702 

! Comment 

here 

0.0000 

! Comment 

here 

1 . 0701 

! Comment 

here 

7.5000 

! Comment 

here 

11 . 9083 

! Comment 

here 

13.3822 

! Comment 

here 

-10.3004 

! Comment 

here 

0.0000 

! Comment 

here 

10.0000 

! Comment 

here 

2 . 8793 

! Comment 

here 

33.8667 

! Comment 

here 

3.2650 

! Comment 

here 

1.0563 

! Comment 

here 

2.0384 

! Comment 

here 

6.1431 

! Comment 

here 

6.9290 

! Comment 

here 

0.0283 

! Comment 

here 

0.0570 

! Comment 

here 

-2 .4837 

! Comment 

here 

5.8374 

! Comment 

here 

10.0000 

! Comment 

here 

1.8820 

! Comment 

here 

-1.2327 

! Comment 

here 

2.1861 

! Comment 

here 

1.5591 

! Comment 

here 

0.0100 

! Comment 

here 

5.1243 

! Comment 

here 

3.0533 

! Comment 

here 

38.9494 

! Comment 

here 

2.1533 

! Comment 

here 

0.5000 

! Comment 

here 

1.0000 

! Comment 

here 

5.0000 

! Comment 

here 
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0.0000  ! Comment  here 
7.7896  ! Comment  here 

7  !Nr  of  atoms;  cov.r;  valency ; a . m; Rvdw; Evdw; gammaEEM; cov . r2 ; 

alf a; gammavdW; valency; Eunder; Eover ; chiEEM; etaEEM; n . u . 


cov  r3;Elp 

; Heat  inc . 

;  n  .  u  . ;  n  . 

u  .  ;  n  .  u  .  ;  n  . 

u . 

ov/ un; vail 

; n . u . ; val3 

, vval4 

c 

1.3742 

4.0000 

12.0000 

1 . 9684 

0 . 1723 

0.8712 

1.2385 

4.0000 

8 .7696 

100.0000 

4.0000 

31.0823 

79.5548 

5.7254 

6.9235 

0.0000 

1.2104 

0.0000 

183.8108 

5.7419 

33.3951 

11 . 9957 

0.8563 

0.0000 

-2.8983 

4 .7820 

1.0564 

4.0000 

2 . 9663 

1 . 6737 

0 . 1421 

14 . 0707 

0.0001 

1 . 9255 

H 

0.6867 

1.0000 

1.0080 

1.3525 

0.0616 

0.8910 

-0.1000 

1.0000 

9.1506 

100.0000 

1.0000 

0.0000 

121 . 1250 

3.8446 

10.0839 

1.0000 

-0.1000 

0.0000 

58.4369 

3.8461 

3.2540 

1.0000 

1.0698 

0.0000 

-15.7683 

2 . 1504 

1.0338 

1.0000 

2.8793 

1.2669 

0.0139 

12 . 4538 

0.0001 

1 .4430 

0 

1.3142 

2.0000 

15.9990 

1 . 9741 

0.0880 

0 . 8712 

1 .1139 

6.0000 

9.9926 

100.0000 

4.0000 

29.5271 

116.0768 

8.5000 

7 . 1412 

2.0000 

0.9909 

14 . 7235 

69.2921 

9.1371 

1 . 6258 

0.1863 

0 . 9745 

0.0000 

-3.5965 

2.5000 

1.0493 

4.0000 

2 . 9225 

1 . 7221 

0.1670 

13.9991 

624.0000 

1.7500 

N 

1.2456 

3.0000 

14 . 0000 

2.0437 

0 . 1035 

0.8712 

1 .1911 

5.0000 

9.8823 

100.0000 

4.0000 

32 .4758 

100.0000 

6.8453 

6.8349 

2.0000 

1.0636 

0.0276 

127 . 9672 

2.2169 

2.8632 

2 .4419 

0.9745 

0.0000 

-4.0959 

2 . 0047 

1.0183 

4.0000 

2 . 8793 

1.5967 

0.1649 

13.9888 

1239.0000 

1.8300 

S 

1 . 9647 

2.0000 

32.0600 

2 . 0783 

0.2176 

1.0336 

1.5386 

6.0000 

9.9676 

5.0812 

4.0000 

35.1648 

112 . 1416 

6.5000 

8.2545 

2.0000 

1 .4703 

9.4922 

70.0338 

8.5146 

28.0801 

8.5010 

0 . 9745 

0.0000 

-10.0773 

2 .7466 

1.0338 

6.2998 

2.8793 

1.8000 

0.0000 

14 . 0000 

180.0000 

2 . 0783 

Si 

2 . 0276 

4.0000 

28.0600 

2.2042 

0 . 1322 

0.8218 

1.5758 

4.0000 

11.9413 

2.0618 

4.0000 

11.8211 

136.4845 

1.8038 

7.3852 

0.0000 

-1.0000 

0.0000 

126.5331 

6.4918 

8.5961 

0.2368 

0.8563 

0.0000 

-3.8112 

3.1873 

1.0338 

4.0000 

2.5791 

0.0000 

0.0000 

0.0000 

180.0000 

2.2042 

X 

-1000.0000 

2.0000 

1.0080 

2.0000 

0.0000 

1.0000 

-0.1000 

6.0000 

10.0000 

2.5000 

4.0000 

0.0000 

0.0000 

8.5000 

1.5000 

0.0000 

-0.1000 

0.0000 

-2.3700 

8 .7410 

13.3640 

0.6690 

0.9745 

0.0000 

-11.0000 

2 .7466 

1.0338 

4.0000 

2.8793 

0.0000 

0.0000 

0.0000 

180.0000 

2.0000 

18 

!  Nr 

of  bonds 

;  Edisl ; LPpen; n . u . 

; pbel ; pbo5 

; 13corr ; 

pbo6 
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pbe2 ; pbo3 ; 

pbo4 ; n . u . 

; pbol ; pbo2 

: ; ovcorr 

1 

1 

141 . 9346 

113.4487 

67 . 6027 

0.1554 

-0.3045 

1.0000 

30.4515 

0.4283 

0.0801 

-0.2113 

8.5395 

1.0000 

-0.0933 

6.6967 

1.0000 

0.0000 

1 

2 

163.6889 

0.0000 

0.0000 

-0 . 4525 

0.0000 

1.0000 

6.0000 

0.5921 

12.1053 

1.0000 

0.0000 

1.0000 

-0.0097 

8.6351 

0.0000 

0.0000 

1 

3 

159.7219 

116.8921 

77 . 9315 

-0 . 4324 

-0 . 1742 

1.0000 

15.0019 

0.5160 

1.2934 

-0.3079 

7 . 0252 

1.0000 

-0 . 1543 

4.5116 

0.0000 

0.0000 

1 

4 

128 . 9104 

171.2945 

100.5836 

-0.1306 

-0 .4948 

1.0000 

26.7458 

0 .4489 

0.3746 

-0.3549 

7.0000 

1.0000 

-0 . 1248 

4 . 9232 

1.0000 

0.0000 

1 

5 

128.7959 

56.4134 

39.0716 

0.0688 

-0.4463 

1.0000 

31.1766 

0.4530 

0 .1955 

-0.3587 

6.2148 

1.0000 

-0 . 0770 

6.6386 

1.0000 

0.0000 

2 

2 

169.8421 

0.0000 

0.0000 

-0.3591 

0.0000 

1.0000 

6.0000 

0 . 7503 

9.3119 

1.0000 

0.0000 

1.0000 

-0.0169 

5.9406 

0.0000 

0.0000 

2 

3 

219.7016 

0.0000 

0.0000 

-0.6643 

0.0000 

1.0000 

6.0000 

0.9854 

5.1146 

1.0000 

0.0000 

1.0000 

-0.0532 

5.1189 

0.0000 

0.0000 

2 

4 

208 . 0443 

0.0000 

0.0000 

-0.3923 

0.0000 

1.0000 

6.0000 

0.3221 

10.5505 

1.0000 

0.0000 

1.0000 

-0.0690 

6.2949 

0.0000 

0.0000 

2 

5 

128.6090 

0.0000 

0.0000 

-0.5555 

0.0000 

1.0000 

6.0000 

0 .4721 

10 . 8735 

1.0000 

0.0000 

1.0000 

-0.0242 

9.1937 

1.0000 

0.0000 

2 

6 

137 . 1002 

0.0000 

0.0000 

-0.1902 

0.0000 

1.0000 

6.0000 

0.4256 

17.7186 

1.0000 

0.0000 

1.0000 

-0 . 0377 

6.4281 

0.0000 

0.0000 

3 

3 

108.9631 

158.3501 

42 . 0558 

0.1226 

-0.1324 

1.0000 

28.5716 

0.2545 

1.0000 

-0.2656 

8.6489 

1.0000 

-0.1000 

6.8482 

1.0000 

0.0000 

3 

4 

85.0402 

118.8680 

75.7263 

0.7080 

-0.1062 

1.0000 

16.6913 

0.2407 

0.3535 

-0.1906 

8.4054 

1.0000 

-0.1154 

5.6575 

1.0000 

0.0000 

3 

5 

0.0000 

0.0000 

0.0000 

0.5563 

-0.4038 

1.0000 

49.5611 

0.6000 

0.4259 

-0 . 4577 

12.7569 

1.0000 

-0.1100 

7 .1145 

1.0000 

0.0000 

3 

6 

191 . 1743 

52.0733 

43.3991 

-0.2584 

-0.3000 

1.0000 

36.0000 

0.8764 

1 . 0248 

-0.3658 

4.2151 

1.0000 

-0.5004 

4.2605 

1.0000 

0.0000 

4 

4 

160 . 6599 

73.3721 

154.2849 

-0.7107 

-0.1462 

1.0000 

12.0000 

0.6826 

0.9330 

-0 . 1434 

10 . 6712 

1.0000 

-0.0890 

4 . 6486 

1.0000 

0.0000 

4 

5 

0.0000 

0.0000 

0.0000 

0.4438 

-0.2034 

1.0000 

40.3399 

0.6000 

0.3296 

-0.3153 

9.1227 

1.0000 

-0.1805 

5.6864 

1.0000 

0.0000 

5 

5 

96.1871 

93.7006 

68.6860 

0 .0955 

-0 .4781 

1.0000 

17 . 8574 

0.6000 

0.2723 

-0.2373 

9.7875 

1.0000 

-0.0950 

6.4757 

1.0000 

0.0000 

6 

6 

109.1904 

70 . 8314 

30.0000 

0.2765 

-0.3000 

1.0000 

16.0000 

0 . 1583 

0.2804 

-0.1994 

8 .1117 

1.0000 

-0 .0675 

8.2993 

0.0000 

0.0000 

10 

!  Nr  of  c 

>ff-diagonal  terms; 

Ediss; Ro; 

gamma; rsigma; rpi; 

rpi2 

1 

2 

0 . 0464 

1.8296 

9.9214 

1.0029 

-1.0000  - 

-1.0000 

0.0000 

1 

3 

0.1028 

1 . 9277 

9.1521 

1.3399 

1 . 1104 

1.1609 

632.0000 

1 

4 

0.2070 

1.7366 

9.5916 

1.2960 

1.2008 

1.1262 

650.0000 

1 

5 

0.1408 

1.8161 

9.9393 

1.7986 

1.3021 

1.4031 

0.0000 
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2 

3 

0. 

0403  1. 

.6913  10. 

4801  0. 

.8774  -1 

.0000  -1. 

0000  0. 

.0000 

2 

4 

0. 

0524  1. 

.7325  10. 

1306  0. 

.9982  -1 

.0000  -1. 

0000  295. 

.0000 

2 

5 

0. 

0895  1. 

.6239  10. 

0104  1. 

.4640  -1 

.0000  -1. 

0000  0. 

.0000 

2 

6 

0. 

0470  1. 

.6738  11. 

6877  1. 

.1931  -1 

.0000  -1. 

0000  0. 

.0000 

3 

4 

0. 

0491  1. 

.7025  10. 

6101  1. 

.3036  1 

.1276  1. 

0173  880. 

.0000 

3 

6 

0. 

1263  1. 

.8163  10. 

6833  1. 

.6266  1 

.2052  -1. 

0000  0. 

.0000 

61 

!  Nr  of  angles ; at 1 ; at2 

; at 3; Thetao, o; ka; : 

kb; pvl ; pv2 

1 

1 

1 

74 . 0317 

32.2712 

0 . 9501 

0.0000 

0 . 1780 

10.5736 

1 . 0400 

1 

1 

2 

70.6558 

14.3658 

5.3224 

0.0000 

0.0058 

0.0000 

1.0400 

1 

1 

3 

65.1700 

8.0170 

7.5000 

0.0000 

0.2028 

10.0000 

1.0400 

1 

1 

4 

65.4228 

43.9870 

1.5602 

0.0000 

0.2000 

10.0000 

1 . 8525 

1 

1 

5 

74 .4180 

33.4273 

1.7018 

0.1463 

0.5000 

0.0000 

1 . 6178 

1 

2 

1 

0.0000 

3.4110 

7.7350 

0.0000 

0.0000 

0.0000 

1.0400 

1 

2 

2 

0.0000 

0.0000 

6.0000 

0.0000 

0.0000 

0.0000 

1.0400 

1 

2 

3 

0.0000 

0.0019 

6.0000 

0.0000 

0.0000 

0.0000 

1 . 0400 

1 

2 

4 

0.0000 

0.0019 

6.0000 

0.0000 

0.0000 

0.0000 

1 . 0400 

1 

2 

5 

0.0000 

0.0019 

6.0000 

0.0000 

0.0000 

0.0000 

1.0400 

1 

3 

1 

72 . 1018 

38 .4720 

1.3926 

0.0000 

0.4785 

0.0000 

1.2984 

1 

3 

2 

89.0416 

36.9460 

0 . 4569 

0.0000 

2.7636 

0.0000 

2 . 0494 

1 

3 

3 

89.9987 

44 . 9806 

0.5818 

0.0000 

0 .7472 

0.0000 

1.2639 

1 

3 

4 

70.3281 

12 . 9371 

7.5000 

0.0000 

0 .7472 

0.0000 

1.2639 

1 

4 

1 

68.3788 

18.3716 

1.8893 

0.0000 

2.4132 

0.0000 

1.3993 

1 

4 

2 

90.0000 

33.6636 

1 . 1051 

0.0000 

0.2638 

0.0000 

1 . 1376 

1 

4 

3 

86.5585 

37 . 6814 

1 .1611 

0.0000 

1 . 7325 

0.0000 

1 . 0440 

1 

4 

4 

74.4818 

12 .0954 

7.5000 

0.0000 

1 . 7325 

0.0000 

1 . 0440 

1 

5 

1 

79.7037 

28.2036 

1.7073 

0.1463 

0.5000 

0.0000 

1 . 6453 

1 

5 

2 

85.9449 

38.3109 

1.2492 

0.0000 

0.5000 

0.0000 

1 . 1000 

1 

5 

5 

85.6645 

40.0000 

2 . 9274 

0.1463 

0.5000 

0.0000 

1.3830 

2 

1 

2 

76.7339 

14 . 4217 

3.3631 

0.0000 

0.0127 

0.0000 

1 . 0400 

2 

1 

3 

56.4426 

17 . 6020 

5.3044 

0.0000 

0.9699 

0.0000 

1 . 1272 

2 

1 

4 

71 . 0777 

9.1462 

3.4142 

0.0000 

0.9110 

0.0000 

1.0400 

2 

1 

5 

63.3289 

29.4225 

2 . 1326 

0.0000 

0.5000 

0.0000 

3.0000 

2 

2 

2 

0.0000 

27 . 9213 

5.8635 

0.0000 

0.0000 

0.0000 

1 . 0400 

2 

2 

3 

0.0000 

0.0019 

6.0000 

0.0000 

0.0000 

0.0000 

1 . 0400 

2 

2 

4 

0.0000 

0.0019 

6.0000 

0.0000 

0.0000 

0.0000 

1.0400 

2 

2 

5 

0.0000 

0.0019 

6.0000 

0.0000 

0.0000 

0.0000 

1 . 0400 

2 

2 

6 

0.0000 

47 . 1300 

6.0000 

0.0000 

1.6371 

0.0000 

1 . 0400 

2 

3 

2 

82.2020 

12 .7165 

3.9296 

0.0000 

0.2765 

0.0000 

1 . 0470 

2 

3 

3 

81 . 1709 

4.2886 

6.5904 

0.0000 

3.0000 

0.0000 

1.2618 

2 

3 

4 

75.9203 

44 . 9675 

0.8889 

0.0000 

3.0000 

0.0000 

1.2618 

2 

3 

6 

83.7634 

5.6693 

2 .7780 

0.0000 

1.6982 

0.0000 

1 . 0400 

2 

4 

2 

55.8679 

14.2331 

2 . 9225 

0.0000 

0.2000 

0.0000 

2 . 9932 
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2 

4 

3 

83.8493 

44 . 9000 

1.3580 

0.0000 

0.5355 

0.0000 

2.5279 

2 

4 

4 

78.7452 

24.2010 

3.7481 

0.0000 

0.5355 

0.0000 

2.5279 

2 

5 

2 

83.8555 

5.1317 

0 . 4377 

0.0000 

0.5000 

0.0000 

3.0000 

2 

5 

5 

97 . 0064 

32 . 1121 

2.0242 

0.0000 

0.5000 

0.0000 

2.8568 

2 

6 

2 

78.3939 

20 . 9772 

0.8630 

0.0000 

2 . 8421 

0.0000 

1.0400 

2 

6 

3 

73.8232 

16.6592 

3.7425 

0.0000 

0.8613 

0.0000 

1.0400 

2 

6 

6 

75.6168 

21.5317 

1 . 0435 

0.0000 

2.5179 

0.0000 

1 . 0400 

3 

1 

3 

71 . 7582 

26.7070 

6.0466 

0.0000 

0.2000 

0.0000 

1 . 8525 

3 

1 

4 

73.7046 

23.8131 

3.9811 

0.0000 

0.2000 

0.0000 

1.8525 

3 

2 

3 

0.0000 

0.0019 

6.0000 

0.0000 

0.0000 

0.0000 

1.0400 

3 

2 

4 

0.0000 

0.0019 

6.0000 

0.0000 

0.0000 

0.0000 

1 . 0400 

3 

2 

6 

0.0000 

31 . 0427 

4.5625 

0.0000 

1 . 6371 

0.0000 

1 . 0400 

3 

3 

3 

84.2807 

24 .1938 

2.1695 

-10.0000 

0 .7472 

0.0000 

1.2639 

3 

3 

4 

84.2585 

44 . 1039 

0.9185 

0.0000 

0 .7472 

0.0000 

1.2639 

3 

3 

6 

73.4663 

25.0761 

0.9143 

0.0000 

2.2466 

0.0000 

1 . 0400 

3 

4 

3 

78.5850 

44.3389 

1.3239 

-26.2246 

1 . 7325 

40.0000 

1 . 0440 

3 

4 

4 

77 . 6245 

32.0866 

1.8889 

-0.9193 

1 . 7325 

0.0000 

1 . 0440 

3 

6 

3 

90.0344 

7.7656 

1 . 7264 

0.0000 

0.7689 

0.0000 

1.0400 

3 

6 

6 

70.3016 

15.4081 

1.3267 

0.0000 

2 . 1459 

0.0000 

1 . 0400 

4 

1 

4 

65.6602 

40.5852 

1 . 8122 

0.0000 

0.2000 

0.0000 

1.8525 

4 

2 

4 

0.0000 

0.0019 

6.0000 

0.0000 

0.0000 

0.0000 

1.0400 

4 

3 

4 

74.2312 

25.7005 

4.3943 

0.0000 

0 .7472 

0.0000 

1.2639 

4 

4 

4 

66.4718 

15.9087 

7.5000 

0.0000 

1 . 7325 

0.0000 

1 . 0440 

6 

2 

6 

0.0000 

31.5209 

6.0000 

0.0000 

1 . 6371 

0.0000 

1.0400 

6 

3 

6 

22 . 1715 

3.6615 

0.3160 

0.0000 

4 .1125 

0.0000 

1.0400 

6 

6 

6 

69.3456 

21.7361 

1.4283 

0.0000 

-0.2101 

0.0000 

1.3241 

31 

u 

o 

-p 

m 

o 

u 

sions ; at 1 ; 

at 2 ;  at 3;  at 4 ;  ;  VI;  V2; 

;  V3;  V2  (BO) 

; vcon j ; n , 

.  u;  n 

0 

1 

1 

0 

0.0930 

18.6070 

-1.3191 

-9.0000 

-1.0000 

0.0000 

0.0000 

0 

1 

2 

0 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0 

1 

3 

0 

1 . 7254 

86.0769 

0.3440 

-4.2330 

-2.0000 

0.0000 

0.0000 

0 

1 

4 

0 

-1.3258 

149.8644 

0.4790 

-7 . 1541 

-2.0000 

0.0000 

0.0000 

0 

1 

5 

0 

4.0885 

78.7058 

0 . 1174 

-2 .1639 

0.0000 

0.0000 

0.0000 

0 

2 

2 

0 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0 

2 

3 

0 

0.0000 

0.1000 

0.0200 

-2.5415 

0.0000 

0.0000 

0.0000 

0 

2 

4 

0 

0.0000 

0.1000 

0.0200 

-2.5415 

0.0000 

0.0000 

0.0000 

0 

2 

5 

0 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0 

2 

6 

0 

0.0000 

0.0000 

0.1200 

-2 .4847 

0.0000 

0.0000 

0.0000 

0 

3 

3 

0 

1.2314 

116.5137 

0.5599 

-4 . 1412 

0.0000 

0.0000 

0.0000 

0 

3 

4 

0 

1.3168 

57 . 0732 

0.2679 

-4 . 1516 

-2.0000 

0.0000 

0.0000 

0 

3 

6 

0 

0.0000 

0.0000 

0.1200 

-2 .4703 

0.0000 

0.0000 

0.0000 

0 

4 

4 

0 

2.0000 

75.3685 

-0 .7852 

-9.0000 

-2.0000 

0.0000 

0.0000 

0 

5 

5 

0 

-0.0170 

-56.0786 

0.6132 

-2.2092 

0.0000 

0.0000 

0.0000 
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0 

6 

6  0  0.0000 

0.0000 

0.1200 

-2 .4426 

0.0000 

0.0000 

0.0000 

1 

1 

1  ] 

0.0000 

48.4194 

0.3163 

-8.6506 

-1 . 7255 

0.0000 

0.0000 

1 

1 

1  2 

>  0.0000 

63.3484 

0.2210 

-8.8401 

-1.8081 

0.0000 

0.0000 

1 

1 

3  : 

5  1.2707 

21.6200 

1.5000 

-9.0000 

-2.0000 

0.0000 

0.0000 

1 

1 

3  L 

l  1.2181 

119.6186 

-1.5000 

-7.0635 

-2.0000 

0.0000 

0.0000 

1 

1 

4  2  -2.0000 

147.2445 

-1.5000 

-7 . 0142 

-2.0000 

0.0000 

0.0000 

1 

1 

4  3  -2.0000 

47 . 8326 

-1.5000 

-9.0000 

-2.0000 

0.0000 

0.0000 

1 

3 

3  ] 

.  -1.8804 

79.9255 

-1.5000 

-4 .1940 

-2.0000 

0.0000 

0.0000 

1 

3 

4  3  2.0000 

96.6281 

-1.5000 

-3.8076 

-2.0000 

0.0000 

0.0000 

1 

4 

4  3  0.1040 

70.1152 

0.5284 

-3.5026 

-2.0000 

0.0000 

0.0000 

2 

1 

1  2 

>  0.0000 

45.2741 

0 .4171 

-6.9800 

-1.2359 

0.0000 

0.0000 

2 

1 

3  l 

l  -2.0000 

156.6604 

1 . 1004 

-7.3729 

-2.0000 

0.0000 

0.0000 

2 

3 

4  3  -0.2997 

152 . 9040 

-1.5000 

-4 . 4564 

-2.0000 

0.0000 

0.0000 

2 

4 

4  3  0.1040 

70.1152 

0.5284 

-3.5026 

-2.0000 

0.0000 

0.0000 

3 

1 

3  3  -2.0000 

22.5092 

1.5000 

-8.9500 

-2.0000 

0.0000 

0.0000 

4 

1 

4  < 

l  -2.0000 

20 . 6655 

-1.5000 

-9.0000 

-2.0000 

0.0000 

0.0000 

9 

!  Nr 

of  hydrogen  bonds ; at 1 ; at 2 ; at 3 

: ; Rhb ; Dehb 

;  vhbl 

3 

2 

3 

2.1845  -2 

.3549  3. 

.0582  19. 

1627 

3 

2 

4 

1.6658  -3 

.8907  3. 

.0582  19. 

1627 

3 

2 

5 

2.6644  -3 

.0000  3. 

.0000  3. 

0000 

4 

2 

3 

1.8738  -3 

.5421  3. 

.0582  19. 

1627 

4 

2 

4 

1.8075  -4 

.1846  3. 

.0582  19. 

1627 

4 

2 

5 

4.0476  -3 

.0000  3. 

.0000  3. 

0000 

5 

2 

3 

2.1126  -4 

.5790  3. 

.0000  3. 

0000 

5 

2 

4 

2.2066  -5 

.7038  3. 

.0000  3. 

0000 

5 

2 

5 

1.9461  -4 

.0000  3. 

.0000  3. 

0000 
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